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Abstract 



The interest in the physical properties of mesoscopic systems experiences a 
gradual increase in recent years. Significant results, which were obtained 
by both theoretical and experimental studies, together with the enormous 
promise of nano-technology applications, contribute to this interest. 

In many mesoscopic systems, the length scale governing the electrons 
motion is small enough to cause quantization of the energy levels. In this 
work we examine two types of such quantum systems: A quantum dot (QD) 
and a nano-wire. In a QD, the system is small enough in all dimensions, so 
that its energy spectrum is entirely discrete, while in a nano-wire the energy 
spectrum is quantized in two of the space coordinates, and continuous in the 
third direction. During the last decade several experimental techniques have 
been developed for manufacturing both kinds of devices, which are currently 
an important tool for understanding low dimensions physics. 

As a result, experiments on mesoscopic samples are conducted nowadays 
in many laboratories, and many interesting results are obtained. For ex- 
ample, some transport properties were recently shown to exhibit interesting 
phenomena in the presence of both interactions and disorder. An analyt- 
ical treatment of these problems is unfortunately difficult, since both the 
disorder and the interactions cannot be considered as a small perturbation. 
A traditional numerical treatment of such systems (exact diagonalization) 
is also problematic, since the dimension of the many-particle Hilbert space 
grows exponentially with the system size. Thus this numerical method is 
usually restricted to small systems. Nevertheless, using analytical methods 
and sophisticated numerical calculations, the interplay between disorder and 
interactions is investigated in the current research for several types of sys- 
tems. 
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ABSTRACT 



In the first half of this work we study a system composed of a QD with 
a single level which is coupled to a one-dimensional (ID) interacting wire 
with spinless electrons. We start by focusing on the fiUing of the dot level, 
the total occupation and the free energy of the system. Using Green func- 
tions technique, we calculate these observables in the non- interacting limit. 
We then investigate numerically, using the density-matrix renormalization- 
group (DMRG) method, two phases of the interacting lead: The Tomonaga- 
Luttinger liquid (TLL) and the charge density wave (CDW) phases. We 
explore the influence of interactions in the lead, as well as dot-lead interac- 
tions, on the width of the dot filling as a function of the chemical potential, 
and on the position of the dot level. In the TLL phase the results are ex- 
plained within the random phase approximation. In the case of a CDW, 
we show that a semi-infinite lead coupled to the dot undergoes a first order 
quantum phase transition when the dot's level crosses the wire's chemical 
potential. 

The Friedel oscillations in the wire, resulting from the dot located at one of 
its edges, are then studied. For the non-interacting case, we develop an exact 
formula of the oscillations in the ID tight-binding model. When interactions 
in the wire are considered, the difference between the two phases of the lead 
are explored. In the TLL phase the oscillations of a clean interacting sample 
decay as a power law, and once a disorder is introduced (Anderson insulator 
(AI) phase), the power law decay is multiplied by an exponential decay term 
due to the disorder. The resulting decay length is shown to increase as a 
function of the interaction strength. On the other hand, when the wire is 
in a CDW phase and the disorder is weak enough, the wire may still be 
described by a Mott insulator phase, and the effect of interactions is the 
opposite. 

We prove that the length scale governing the exponential decay, in the 
AI phase, may be associated with the Anderson localization length and thus 
be used as a convenient way to determine the dependence of the localization 
length on disorder and interactions. Our results show a decrease of the 
localization length as a function of the interaction strength, in accordance 
with previous predictions. 

In the second half of the research we study two cases of an isolated two- 
dimensional QD. We begin by exploring some properties of a disordered QD 
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consisting of interacting spinless electrons. Since the size of the relevant 
Hilbert space is huge, such a problem cannot be solved by an exact diag- 
onalization (except for very small systems). We thus use the sophisticated 
particle-hole DMRG (PH-DMRG) method, showing that its approximation 
for the ground-state energy is much more accurate than that of the Hartree- 
Fock (HF) method. We also suggest an improvement of the PH-DMRG 
truncation algorithm, which reduces the error rate of the traditional method 
by almost 30 percents. 

As an apphcation of the improved PH-DMRG method we calculate the 
addition spectrum of the QD. We present the improvement of the PH-DMRG 
results comparing to the HF approximation in three aspects: the error rate, 
the average and the fluctuations of the addition spectrum. 

Finally we study the magnetization of a QD with spin 1/2 electrons, in 
the presence of spin-orbit (SO) coupling. We calculate the g-factor and the 
expectation values of the spin operators in the ground state. We find that 
when the dot is occupied by an even number of electrons, there is a level 
crossing between the two lowest many-body eigenfunctions as a function of 
the SO scattering rate, resulting in a finite magnetization of the ground 
state. This is a clear signature of the interplay between SO scattering and 
interactions, and may have a significant influence on g-factor measurements. 
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Chapter 1 
Introduction 



In recent years there is a gradually increasing interest in the physical proper- 
ties of mesoscopic systems. Significant results, which were obtained by both 
theoretical and experimental studies, together with the enormous promise of 
nano-technology applications, contribute to this interest. 

Mesoscopic physics deals with systems whose sizes are between the mi- 
croscopic and the macroscopic regimes [U 12] . Such an intermediate regime 
has a length-scale which is small enough so that quantum physics phenomena 
are relevant, and still large enough to enable laboratory experiments. This 
unique combination, gives a challenging topic for research, and attracts many 
recent investigations. 

The classification of a physical system as a mesoscopic one is usually done 
by an examination of the electron motion, checking whether it is coherent, 
i.e., if the electron's phase is conserved |3|. The loss of phase conservation can 
be caused by inelastic collisions with phonons or with other electrons, thus 
one can say that for mesoscopic systems, the mean distance which the elec- 
tron passes between such subsequent inelastic collisions, is significantly larger 
than the system size. Since the probability for inelastic collisions decreases 
when the temperature is lowered, the recent advances in low-temperature 
physics make experimental development of mesoscopic systems feasible. 

It is important to notice, however, that there is another important type 
of collisions, which does not cause a phase decoherence. For example, when 
an electron experiences an elastic collision with an impurity, its phase is 
conserved, so that the motion is coherent. Such collisions, therefore, can be 
used to further classify the system. Usually one defines the electron's mean 
free path as the mean distance which the electron passes between subsequent 
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elastic collisions with impurities. Denoting the mean free path by /, the ratio 
between / and either the system size L or the inverse Fermi momentum kp^ 
can characterize the motion. For / > L the electron is hardly influenced by 
the impurities, so that its motion is ballistic. In the other limit, when there 
are lots of impurities, i.e., the disorder is strong, so that / < kp^ (where 
kp is the Fermi momentum), the electron's motion is limited in space, and 
the electronic wave functions are localized. In the intermediate regime, for 
kp^ < I < L, the motion is diffusive [3]. 

In many mesoscopic systems, the length of one (or more) of the dimensions 
is small enough to cause quantization of the energy levels. In our research 
we will examine a quantum dot, in which the system is small enough in all 
dimensions, so that its energy spectrum is entirely discrete. In some parts 
of this work, we will be interested also in a one-dimensional lead which is 
connected to the dot. Such a device, composed of a dot and one or two leads, 
can be used for various transport and charge sensing measurements, and can 
be an important tool for understanding low dimensions physics. 

In this chapter we shortly describe the two main parts of the physical 
model. We first briefiy introduce the quantum dot, after which we will take 
a short tour in the one- dimensional world. We conclude this chapter with 
the motivation for the research, and the outline of the rest of the thesis. 

1.1 Quantum Dots 

A quantum dot (QD) is a mesoscopic physical device, in which the motion of 
the electrons is restricted in space, so that the energy levels are quantized in 
all three dimensions. As a result, the QD has a completely discrete energy 
spectrum [H O [6l [7j . 

Experimentally, there are several physical realizations, such as metallic 
or semiconducting nano-particles, in which the electrons are restricted to a 
certain region in space. Most experimental QDs are based, however, on semi- 
conductor heterostructures. A two-dimensional (2D) plane, into which the 
electrons are confined, is created, and by further applying a voltage above 
or below specific points in the plane, one is able to restrict the electrons' 
motion inside a small region in this 2D plane. An example for such a con- 
struction is shown in Fig. II. 1[ Therefore, a common QD situation is that of 
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electrons moving in a confined 2D area, which is separated from other regions 
by potential barriers. 

Figure 1.1: A typical example of a 2D QD (drawn schematically on 
the left, while the actual electron micrograph is shown on the right). The 
electrons are confined vertically to the 2D interface between GaAs and 
AlGaAs, and negative voltages applied to the metallic gates confine them 
also laterally. The transport of electrons through the dot is indicated by 
red arrows. The mean free path and the coherence length are of the order of 
10 microns, so that the electronic motion is coherent and ballistic (picture 
taken from Ref. |7j). 

1.1.1 Hamiltonian for a QD 

As many QDs are defined on a 2D electron gas, one of the convenient ways 
to model them numerically is to use the Anderson model, which is derived 
from a tight-binding description of a 2D dot. Assuming that an electron is 
restricted to lattice sites, with the ability to hop from one site to another 
(usually limited to be one of its neighbors), one can therefore replace the 
continuous space with an effective 2D lattice model, for example of A rows 
and B columns, which sum up to AB sites. 

For each one of the AB sites there is a well defined on-site energy, e, which 
is the potential that the electron feels when it is located on that site. In order 
to model impurities in the QD, a common method is to take a different on- 
site energy for each site. These energies are usually determined randomly, 
using a uniform distribution between [—W/2,W/2], where W describes the 
strength of the disorder. Therefore, by using the fermionic operators a] and 
di for a creation and annihilation, respectively, of an electron at site i, the 
on-site energy term in the Hamiltonian can be written as Yli^i^l^i- 

The hopping element will be written, in a similar notation, as —t ^-^ (a|dj+ 
H.C.), where denotes nearest neighbor (NN) sites i and j, and the over- 
lap integral between NNs is represented by t and is taken to be constant. 
Finally, we get the following Hamiltonian: 
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Spin-Orbit Coupling 

For spin 1/2 particles the Hamiltonian Eq. (11. ip should be slightly modified 
in order to describe the physics of the system. It is convenient to separate the 
site index from the spin index, and to use two indices: i which describes the 
site index, and a = i which represents the spin projection. For a system 
of sites with spin 1/2 electrons, one should thus diagonalize a matrix of 
size 2N X 2N. 

In some cases there is no spin-orbit coupling, i.e., the orbital and the spin 
degrees of freedom are completely decoupled. In these cases, the system, and 
the resulted Hamiltonian matrix, has an orthogonal symmetry. Effectively, 
the Hamiltonian remains identical to Hqj^ described above, and with the new 
indices notation it can be written as 



However, when spin-orbit coupling is taken into account, the system's 
symmetry is no longer orthogonal, but rather symplectic ^ t9j . Such a cou- 
pling can result, e.g. from the magnetic moments of the impurities, allowing 
an electron to change its spin projection during hopping between sites. The 
Hamiltonian for such a system can be thought of as an x matrix, in 
which every element is a quaternion with 4 degrees of freedom. Diagonaliza- 
tion of such a Hamiltonian results in pairs of degenerate eigenvalues, which 
is known as the Kramers degeneracy |10j . 

In practice, in common QD devices (see e.g. Fig. II. ip . the electrons 
are confined to a 2D interface separating between regions having different 
electric potentials. As a result, a perpendicular electric field is affecting the 
electrons, leading to the spin-orbit Rashba term [11] in the Hamiltonian, 
HRashba = aE^ia x p)^. 

A useful method for considering such a term for a lattice, which was 
suggested by Ando [12] , is by replacing the overlap integral t of Eq. (I1.2p by 
two 2x2 matrices \4 and Vy which encompass the probabilities to change or 
preserve the spin projection in transitions in the x and y directions. If we now 
denote the 2D-lattice site with row and column indices m, n, respectively, the 
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resulting Hamiltonian reads 

"^QD ~ ^m,nO'm,n,a^rn,n,a (1-3) 

m,n,a 

m,n,a,cr' 

The matrices written by Ando are composed of the identity matrix and 
Pauh matrices (Jx,cry with the prefactors Vi, for a motion without changing 
the spin, and V2, for spin-flips. Note that the Rashba term leads to an 
influence of on the motion in the y direction, and of ay in the x direction. 
The matrices are thus 

v^-(l !n : y,-(^i (1.4) 



A correct building of the V matrices is essential in order to conserve 
the system symplectic symmetry. The main manifestation of such a system 
symmetry is a sign flip under a 27t rotation, preserving a regular symmetry 
under rotations of Air. 

Magnetic Field 

In order to insert an in-plane magnetic field dependence into the Hamiltonian, 
one adds the term Hb = <? ■ (fisH) to the Hamiltonian, where a represents 
a vector of Pauli spin matrices, and Hb = 2^ Bohr magneton. In 

the presence of spin-orbit coupling, the in-plane magnetic field yields a non- 
diagonal term. For example, for a magnetic field in the x direction, i.e., 
H = Hx, one has Hb = ^bHcTx- As a result, for the i-th energy level, the 
magnetic field couples the elements |^, t> and |z, i> through 

A Perpendicular magnetic field gives a simpler term, since is diagonal 
in the spin indices t and |. However, in that case one should consider also 
orbital effects and, as a result, the hopping term we take includes the dia- 
magnetic coupling through the Peierls substitution p!3]. Taking a Landau 
gauge for the vector potential A = Hyx one gets a revised hopping element, 
since in the x direction [T3j the hopping element t, or the matrix V^, if spin- 
orbit coupling is considered [151 HB] , should be multiplied by exp ( i 



<t>0 

where m is the row number, s is the lattice constant, and 0o = hc/e is the 
magnetic flux unit. 



6 



CHAPTER 1. INTRODUCTION 



Interactions 

When the interactions between electrons cannot be neglected, an appropriate 
term should be added to the Hamiltonian. Since in metallic materials there 
is a substantial screening effect, one usually restricts the interaction terms to 
have a short range. For spinless electrons it is thus sufficient in many cases 
to deal with NN interactions, so that the interaction term can be written as 

^iiT^ = f/„n ^ aj a] ttj Oi , (1.5) 

while in other cases one may consider the full long-ranged Coulombic term, 
i.e., 

^{Coulomb) = u^Y^ ^ala]ajai, (1.6) 

where r^j, the distance between the sites i and j, is measured in units of the 
lattice constant. 

When the spin of the electron is also considered, the most important 
interaction term is the on-site interaction between spin-up and spin-down 
electrons. This is the Hubbard term which is written as 

^{Hubbard) = u^Y^ al^di^^al^tti,^. (1.7) 

i 

The innocent look of the interaction terms hides the huge impact they 
have on the calculation of the system's physical properties. This results 
from the fact that the system states are multi-particle states, and thus the 
dimension of the Hilbert space grows exponentially with the system size. 
This issue is described in detail in chapter [21 

1.1.2 Coulomb Blockade 

Once the QD is connected to two metallic leads, applying a bias voltage 
on the leads can cause an electric current through the dot. An interesting 
phenomena occurs when the QD is weakly coupled to the leads, so that 
the transport of electrons to and from the dot is by tunneling. Usually, 
tunneling of an electron into the dot is blocked by the Coulomb repulsion 
of the electrons which are already inside the dot. However, by changing 
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the gate voltage one can compensate for that energy difference, and for an 
appropriate value of Vg the number of electrons in the dot can increase by 
one, and a peak in the conductance through the dot appears. 

This phenomenon, named the Coulomb blockade, was first observed [17] 
as early as in 1968. The gradual increase of the gate voltage was found to 
cause jumps in the current through the quantum dot in specific values of the 
gate potential. These jumps are easily seen in I-V (or vs. V) curves. The 
energy values in which the jumps occur were shown to be sample-dependent, 
but subsequent measurements of the I-V characteristics for the same sample 
gives the same values each time. 

For zero temperature, the conductance is possible only through one of 
the empty discrete energy levels in the dot, e^. The transition is allowed 
only if the total energies when there are — 1 or electrons in the dot are 
identical. The electrostatic energy in the dot can be written in the constant 
interaction (CI) model [3] as U{N) = -cNVg + e^N{N - 1) /2C, where C 
is the total capacitance between the dot and the leads. Let's suppose that 
there are currently — 1 electrons in a dot. The condition for the transition 
of the iV-th electron into the dot is 

N-l N 

J2e, + U{N-l) = J2(^i + U{N). (1.8) 

1=1 1=1 

Substituting the formula for U (N) given above, leads to the relation which 
determines fiN = eVg, the chemical potential of the dot in which the conduc- 
tance peak occurs, as 

= eVg = eN + {N- l)eVa (1.9) 

After the A^-th electron was added to the dot, a further increase of Vg 
blocks again the transport channel, until the next electron will be in the 
appropriate position to enter the dot. As a result, peaks will appear in the 
curve of as a function of the gate voltage, for specific values of Vg. In 
practice, one can use this method in order to obtain a control on the number 
of electrons in the dot. Starting from an initial state in which the dot is not 
charged, one can gradually increase the gate voltage and count the current 
peaks, and thus know exactly how many electrons have transferred into the 
dot. 
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1.1.3 Addition Spectrum 

Using Eq. ( II. 9p it is easy to define the addition spectrum as the change in 
chemical potential required to add the A^-th electron to the dot, i.e., the 
distance between the subsequent peaks — 1 and A^: 

Ag^^ = fiN - f^N~i = Cat - eN-1 + e^/C. (1.10) 

Thus, based on the assumptions of the CI model, one finds that the addi- 
tion spectrum consists of the level spacing with an extra charging energy. The 
reason that the distance between peaks seems sometimes constant, is that in 
most cases the charging energy is very large comparing to the level spacing 
fluctuations. However, by subtracting this constant charging energy from the 
measured distances, one can reconstruct the probability distribution of the 
level spacings. For disordered quantum dots in the diffusive regime, the level 
spacings have a well known probability distribution, the Wigner (or Wigner- 
Dyson) distribution, which was first predicted in the 1950's in calculations 
of energy levels in the nucleus [HI [8]. The distribution has a different form 
for different symmetry classes (orthogonal, unitary and symplectic), which 
are sometimes denoted by GOE, GUE and GSE (for Gaussian orthogonal 
ensemble, and similarly for the other two). The Wigner-Dyson distribution 
has become an important part of the random matrix theory [9], [19] . 

Unfortunately, several experiments on QDs which were performed in the 
last decade are not consistent with the results of the random matrix theory 
[201 [211 [22] . Foi' example, while it was found that the mean level spacing can 
be described by the CI model, the A2 distribution does not fit the predicted 
Wigner-Dyson distribution. It is thus clear that for obtaining an accurate 
description of some measurements, the effect of interactions should be con- 
sidered beyond the CI model. 

In order to take into account the interactions in the system, one needs to 
define the addition spectrum in an alternative way, without any assumption 
on the interactions. In general, the gate energy in which the iV'th conduc- 
tance peak occurs must be equal to the difference in the ground-state energies 
for — 1 and A^ electrons. Therefore, fiN = Egs{N) — Egs{N — 1), which 
leads to the exact definition of the addition spectrum as 

Af ) = /i^ - = ^^.(A^) - 2EgsiN - 1) + ^^.(A^ - 2). (1.11) 



1 . 2. ONE-DIMENSIONAL LEAD 



9 



Such an expression requires the knowledge of three different ground-state 
energies, with consecutive electronic occupation. 

1.2 One-Dimensional Lead 

The previous section mentioned some of the physical results of connecting the 
QD to external leads: Current can now be transferred through the QD; The 
coupling of the dot levels with the leads causes a broadening of the energy- 
levels width; If this coupling is weak enough the discrete levels of the QD can 
still be seen experimentally. Of course, more complicated properties of the 
QD can be measured as well. However, it is interesting enough to explore 
these simple phenomena for different kinds of the external leads. What will 
be the influence of the lead on the QD when the lead is not described by 
a simple Fermi liquid theory? Does it have a similar influence, or, perhaps, 
new physics may be explored? 

For answering these questions, we would like to connect a special kind of 
external lead to the QD: a one-dimensional (ID) lead. A device is defined 
as ID if its energy spectrum is quantized in two of the space coordinates, 
while it is continuous in the third direction. During the last decade several 
experimental techniques have been developed for manufacturing various ID 
devices, such as carbon nanotubes [231 EH ESI ES, EZ] , polymer nanofibers 
[28] and semiconducting nanowires [29l [30]. For example. Fig. 11.21 shows a 
silicon quantum wire, which is suspended in a highly doped silicon film in a 
silicon-on-insulator substrate. 

Figure 1.2: SEM micrograph of a suspended silicon quantum wire in a 
highly n-doped silicon-on-insulator film. The wire width is 80nm and the 
length is 1.5fim (taken from Ref. [29j). 

It is thus clear that the creation of a ID lead is nowadays experimentally 
possible. As a result, a door to a new world of physical models is opened. In 
this section we give a brief introduction to that world. 

1.2.1 Hamiltonian for a ID Lead 

We start by presenting the Hamiltonian of a ID lead. As in the QD case, we 
restrict ourselves to the tight-binding description of the electronic orbitals. 
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and thus we start the discussion from a ID lattice. The electrons' motion is 
thus not continuous, as they jump from one site to another. The potential 
energy of the electrons originates from the on-site energy of each lattice 
site. If every site has the same on-site energy, all energies can be rescaled 
with respect to that energy, and it can be excluded from the Hamiltonian. 
However, when there is a difference between lattice sites, e.g. when the unit 
cell consists of two or more atoms, or when impurities are involved, this term 
must be included. 

The kinetic energy of the electrons is written through the hopping el- 
ement, which is in charge for the electrons' motion along the lead. As in 
the QD case, we restrict this term to the physical intuition of NN hopping 
only. This, of course, has a general justification since the hopping probability 
decays exponentially with distance. 

As a result, for electrons with spin but without spin-orbit couplinjl], one 
gets the following Hamiltonian: 

Aead = - i + H.C.), (1.12) 

where c]^ (q.o-) is the creation (annihilation) operatoj^ for an electron with 
spin a in lattice site number i, which has an on-site energy ej. For spinless 
fermions the Hamiltonian is even simpler, and has the form 

^lead = E '^^^^ - ^ + ^•^•)- (1-13) 

i i 

The addition of interactions is straightforward. As noted in the previous 
subsection, usually one can consider only the most important term, which 
is the Hubbard interaction in the spin 1/2 case, and the NN interactions 
for spinless fermions. From now on we will focus on the spinless fermionic 
case, without the existence of a disorder. Therefore, if we assume a positive 
background in the lattice sites, one can write the interaction term as 

^^ad = IY.(4c^ - 1/2)(c,Vq+i - 1/2), (1.14) 

i 

^Spin orbit coupling in the lead is not discussed in this thesis. 

^For convenience, along this thesis the dot operators are denoted by d whereas the lead 
operators are denoted by c. 



1 . 2. ONE-DIMENSIONAL LEAD 



11 



and combining it with the previous term results in the form 

^lead = -t5^(cjQ+i +ff.C.) +/5^(c1q - ^){cl+A+l " (1-15) 
i i 

Note that since multiplying the entire Hamiltonian by a constant factor 
does not change the physics of the problem, practically this Hamiltonian 
has only one free parameter, which is the ratio I/t (given that t is taken as 
positive). 

1.2.2 ID Spin Chain 

Interestingly, the spinless fermionic interacting Hamiltonian which is de- 
scribed by Eq. fll.lSp . can be shown to be equivalent to an entirely different 
type of ID system, that of a spin chain [311 [32]. Let's assume we have a chain 
of spin 1/2, i.e., we have a chain in which there is a spin Si on every site. The 
spin Si is defined by the relation Si = hai/2, where (Xj are the Pauli matrices, 
and the three components of the spin have the regular angular momentum 
commutation relations. 

The Heisenberg Hamiltonian for such a chain considers only NN inter- 
actions between spins, with a rotational symmetry in all directions. The 
Heisenberg Hamiltonian is written as 

-^Heisenberg = ^ ] ^*S'jS'j-|_i = ^ ^ ^i^i ^i+l ~^ Sf Sf^^ + S^ Sj^^^). (1.16) 

i i 

When J is positive the resulted ground state is anti-ferromagnetic, whereas 
for a negative sign of J, a ferromagnetic state is favored. 

When the rotational symmetry is broken, however, the Hamiltonian is 
modified to have different coupling constants (J's) in each direction. A fa- 
mous case involves symmetry breaking in the z direction, whereas the system 
is still invariant for rotations in the x - y plane. The result, known as the 
XXZ Hamiltonian, is thus 

Hxxz = ^ ^Jxy{SfSi_^^ + S'^Sf_^_^) + JzS-S-^^ , (1-17) 

i 

with A = Jz/Jxy as a single free parameter . 

We now wish to show the equivalence of Eqs. (Il.isp and (11.17^ The 
first step is to relate the spin-down state to an empty site in the fermionic 



^In the next steps we follow the route of Ref. [3T] 



12 



CHAPTER 1. INTRODUCTION 



language, and the spin-up state with an occupied one. The dimension of the 
Hilbert space is exactly the same since each site contributes 2 possibilities 
in both representations. Such a mapping can be easily obtained by defining 
— c\ and = c\ci — 1/2, and it is easy to show that the commutation 
relations in each site are obeyed. However, the inter-site commutation rela- 
tions of the spin chain and the anti-commutation relations of the fermionic 
case are different. 

In order to preserve the anti-commutation relations without destroying 
the local commutation relations, Jordan and Wigner has proposed the map- 
ping 

St = 4 exp(z7r J] cjc,), (1.18) 

St = 4q-1/2. 

Using the Jordan- Wigner transformation, one can rewrite the terms of 
the XXZ Hamiltonian in the fermionic language. For example, the x - y 
term in the Hamiltonian can be written as 



Qx cx j_ dy dy _ cx cx , dy cv n iQ^i 



for which it is easy to substitute the Jordan- Wigner mapping, and to get 

Sti^r = cj+iexp(i7r J2 exp(-i7r ^c]c^)q (1.20) 

j<i+l j<i 

= cj+i exp(i7r4cj)cj 

= cl^i{l + inclci + . . .)ci 

since all the other elements end with CjCj and thus vanish. 
In a similar way one gets S^^S^ — clci^i, and thus 



H 



XXZ 



-Y^4+iCi + H.c.) + j,{cl^^Ci+i - ^)(cJq - ^; 



(1.21) 



The final step is to make a canonical transformation in which the mo- 
mentum of the fermions is shifted by tt, i.e., we multiply Cj by a factor 
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exp(i7rj) = (— l)-'. With such a transformation the first term in Eq. f ll.2ip 
gets a minus sign, so finally 



Hxxz — 



Jxy /+ lx/+ 1 

y 



1.22) 



Taking t = ^ and / = one gets to the original fermionic Hamiltonian 
Eq. (11.151) . We thus conclude that the XXZ spin chain and the spinless 
fermionic system have identical physical properties. The significance of that 
conclusion will become clear shortly. 



1.2.3 Phase Diagram 

As noted in the case of the QD Hamiltonian, when interactions are included 
in the Hamiltonian, the system in general cannot be exactly solved. How- 
ever, for some special cases (many of them low- dimensional systems) specific 
methods were developed, so that their physical properties can be found ex- 
actly. 

One such method for exactly solving ID systems is the Bethe ansatz 
technique. It was first suggested [33] by Bethe, in 1931, for the XXZ spin 
chain problem with periodic boundary conditions [M], and since then it was 
extended and used in more complicated systems. As a result, the physical 
properties of the XXZ spin chain, as a function of A = Jz/Jxy, are known 
exactly [35] , and we now review the main results which are important for our 
purpose. 

We first point out that the transformation J^y — Jxy and Jz Jz of 
the XXZ spin chain is identical to replacing the operators , Sf by (— l)*^^^' 
and {—lySf. Therefore, when taking A — —A, the physical system remains 
almost the same, except for a change of the spin-spin coupling between anti- 
ferromagnetic and ferromagnetic coupling^]. 

When A is exactly —1 or +1 the system is isotropic, whether in a fer- 
romagnetic or an anti-ferromagnetic state. Different phases of the system 
evolve for the three regimes divided by these points: A<— 1,— 1<A<1 
and 1 < A. 

■*Note that the transformation J^y — > — Jxy is accompanied, in the fermionic language, 
by an inversion of the energy band, because t changes its sign. 
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For A < — 1 the XXZ spin chain is in a ferromagnetic phase. In the 
ground state of such a system all the spins point to the same direction, 
usually defined along the z axis. Obviously, without a magnetic field, the 
ground state is doubly degenerate with Sz = ±N/2. The first excited states, 
which can be exactly found by the Bethe ansatz method, are magnons with 
Sz = ±{N/2 — 1), and there is an energy gap which is linear in |A| — 1. In 
the limit A — ^ — 1, therefore, the system becomes gapless. 

The 1 < A phase is similar to the A < — 1 one, but instead of a ferromag- 
netic state the ground state is anti- ferromagnetic (Neel). The ground state 
is doubly degenerate regarding the two possible anti-ferromagnetic orders, 
TiTi ■■■ itit and the excited states produce a gap which grows as 
|A| — 1. Again, in the limit A ^ 1 the system becomes gapless. 

The interesting regime is the intermediate one, — 1 < A < 1, for which 
a full solution using the Bethe ansatz technique shows that the system is 
gapless. This phase is called the XY phase. We remark that this phase 
includes the non-interacting point A = 0. The three different phases are 
schematically shown in Fig. 11.31 



- 1 1 A = J /J 

Ferromagnetic . XY , AFM ^ 
^ 1 > 



Phase Separation I XLL ' CDW 

-2 2 



I/t 



Figure 1.3: The phase diagram of the spin chain and the spinless 
fermionic lattice. The different phases are shown as a function of the 
parameter A = Jz/Jxy in the spin chain case, and I/t in the fermionic 
case. 



We now want to understand these results in the fermionic case. We use 
the equivalence shown in the previous section, in which the spin-up (down) 
case in the spin chain was shown to be identical to an occupation (absence) of 
a fermion. The XXZ ferromagnetic phase tells us that the fermionic ground 
state for I < —2t is either entirely occupied or entirely empt}0. The energy 
gap points out that this phase is an insulator. 

^When a system is explicitly enforced to have a certain filling factor (e.g. to be half 
filled) the electrons prefer to gather, and the result is called a "phase separation". 
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Similarly, the XXZ anti-ferromagnetic phase is also related to an insulat- 
ing fermionic phase occurring for / > 2t, which is the charge density wave. 
The degenerate XXZ ground state is related to the two possible charge den- 
sity wave orders in which alternating sites are occupied, either like • o • o ■ • • 
or o • o • ■ ■ ■ . 

The middle phase, the gapless XY phase in the XXZ chain, is projected 
onto a metallic fermionic phase, with —2t < I <2t. It can be shown that the 
XY phase can not be described by the Fermi liquid theory; Rather it belongs 
to a different universality class, called Tomonaga-Luttinger liquid. This class 
will be briefly introduced in the next section. 

1.3 A Short Introduction to Tomonaga-Luttinger 
Liquids 

The inclusion of electron-electron interactions in the Hamiltonian has usu- 
ally a very significant effect on the possibility to exactly calculate physical 
quantities of the system. These difficulties are in general both theoretical 
and numerical. However, a physical interacting system may sometimes be 
equivalent to a non-interacting system, and thus can be exactly analyzed. 
A good example for that is the case of interacting electrons in two or three 
dimensions, a system which is usually referred to as a Fermi liquid. Landau, 
in his famous Fermi liquid theory [3BJ, has proven that instead of calculating 
directly the Fermi liquid properties, one can use its equivalence to a Fermi gas 
of quasi-particles for which these properties can be easily found [H7]. These 
quasi-particles can be shown to share the same important physical proper- 
ties, such as charge and spin, with the original electrons, so that instead of 
the difficult task of solving a system of interacting electrons, one can more 
easily solve a system of quasi-electrons gas. 

Unfortunately, in one dimension it can be shown [381139] that an excitation 
which is composed of quasi-particles is not stable, so that quasi-particles 
cannot describe the ID system. The stable excitations, on the other hand, 
do not have the physical properties of the electrons. In general, the Fermi 
liquid theory cannot be used in one dimension. 

A significant progress in the ID world was obtained by the works of 
Tomonaga [10] and Luttinger [H], and we now briefly review their main 



16 



CHAPTER 1. INTRODUCTION 



steps. The first step is the linearization of the dispersion relation of a free 
electron, which moves in a ID system of length L with periodic boundary 
conditions. It is clear that the energy levels which contribute to the transport 
properties are those which are close to kp, so that one usually neglects the 
levels which are far below or far above kp. Since in one dimension a Fermi 
surface does not exist, but there are only two Fermi points {^kp) instead, the 
contributing electrons are split to "right-movers" {k ~ kp) and "left-movers" 
{k ~ —kp), and thus one can write up to first order in k — \kp\ 



de 

+ -j^\k=+kF\k — kp) 
ep + —\k=-kp{k + kp) 



[1.23) 




Figure 1.4: The dispersion relation e ~ /c^ and the approximation by the 
right and left branches near zizkp. The states taken by Tomonaga's model 
are schematically signed by the red crosses, while the states added by the 
Luttinger model are colored blue. 



These two branches are the linear lines shown in Fig. 11.41 which can be 
used to approximate the quadratic dispersion relation in the vicinity of the 
Fermi points. Since jT\k=±kF = ^'Vp, one gets 



e± = ep + Vp{±k — kp). 



[1.24) 



This energy dependence can be denoted by the following fermionic Hamil- 
tonian, which neglects the constant Fermi energy. 



H = ^ vp{rk - kp)dl f^dr,k, 



[1.25) 



k,r=±l 
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where r = ±1 is used to discriminate between right and left movers, and alj, 
{dr^k) denotes the creation (annihilation) of an electron occupying the k-th 
energy level with a direction r. 

In order to proceed from Eq. (11.251) it is clear that the first question is 
what are the values that k can take. If we restrict k to be in the vicinity 
of ±kp, where the approximation is rigorously valid, the model cannot be 
exactly solved. Tomonaga was the first to offer a creative idea: Since the 
states which are far below kp should not change the basic physics of the 
system, one can include these states in the model, so that it will become 
easier to solve. He offered to include, in the sum over k, all the states which 
are either on the right-moving branch with > 0, or on the left-moving one 
with A; < 0. Such states are schematically drawn in red in Fig. II. 4[ While 
this slight change in the model is still not enough in order to exactly solve 
the Hamiltonian, yet an approximate solution can be found, as was shown 
by Tomonaga himself. 

The next significant step was suggested by Luttinger. His idea was to 
include the bottom of the two branches entirely, i.e., take the right branch 
for — oo < k < kp, and the left branch for —kp < k < oo (i.e., add also 
the blue states in Fig. II. 4p . With this change the model becomes exactly 
solvable, and nowadays the model is mainly named after Luttinger. 

We now present the main steps required to solve the Hamiltonian. First 
we define an operator (for each branch), which changes the electrons momen- 
tum by k: = J2q^l+k^g- turns out that in the right branch with 
A; > creates a particle-hole pair, while a k < term annihilates it, and the 
opposite happens in the left branch. Since pk is built out of two fermionic 
operators, it is not surprising to find out that its commutation relations are 
bosonic, up to a multiplicative factor. It is thus convenient to exactly define 
the bosonic operators bI and Bk as a function of p±k, so that [Bk, bI,] = 6k,k'- 

Defining two number operators, which count the particles on the right 
and left branches, as N^. = po = Ylin'^ln'^r,ni and using the Kronig identity, 
one can express the Hamiltonian Eq. (I1.25p . as a function of the bosonic 
operators, as 




(1.26) 
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This Hamiltonian represents three types of excitations. The first term 
is simply a momentum change of an electron in one of the branches. The 
second term results from two other optional excitations: the jump of an 
electron from one branch to the other (current excitations), and the addition 
of an even number N of electrons to the system, N/2 to every branch (charge 
excitations). 

By transforming the density operators to real space, one can get a more 
convenient representation of the Hamiltonian. We write 

P^(^) = $^exp(-ifca:)p± = iV± + ^ exp(-i/ca:)p± (1.27) 

k fc^O 

and with these definitions and a bit of simple algebra one can show that 

H dx[p+{x)p+{x)+ p'{x)p-{x)]. (1.28) 

Looking at the form of Eq. fll.28p it is obvious that the Hamiltonian can 
be easily written as a function of the "charge density" p^{x) + p~(x), and 
the "current density" p^{x) — p~{x). This can be obtained by definingj the 
fields </)(x) and 6{x): 

0(x) = -[N^ + N_)—--Y, \ {Pt+Pp) (1-29) 

e{x) = {N^ + N^)— + -Y,-\^{pt-Pv)- 

With this definition one finds out that \/(f){x) = — ^[p+(x) + p^{x)] and 
\/9{x) = jIp'^^x) — p~(x)], leading to the nice form 

H=^[ dxvF[{s^(P{x)f + {s^e{x)Y]. (1.30) 

^TT Jo 

The advantage of this Hamiltonian form, is that it is a free-particle 
bosonic Hamiltonian so that its solution is exactly known. However, since we 
haven't yet considered the electron-electron interactions, it may not be such a 
surprise. But, as we now show, the great importance of the Luttinger model 
is that up to some prefactors, interactions do not change this Hamiltonian 
form. 



^The two fields ip{x) and 0{x) are named and defined in different ways by different 
authors. We adopt the version of Ref. [3T]. 
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For spinless electrons there are two kinds of interactions, which are his- 
torically denoted as g2 and g^. In Fig 11.51 these processes are showr]l|. 




Figure 1.5: The interaction process g2 and 34 which are considered in 
the Luttinger model. The first describes a momentum transfer between 
electrons from different branches, while the second describes a transfer of 
momentum between electrons in the same branch (different branches are 
presented in different colors). 



The Hamiltonian terms of these interaction processes are easily written 
in the same form of the non-interacting Hamiltonian: 



1 94 \^ \^ t t ^94 

> > a , ,ar,a ,aa = 

k p,q 



= ijJ^Yl = 1^; J^'^ptp-k + Pk p~k)^ (1-31) 

k p,q k 



k p,q k 

After the transformation to real space one gets 

K = J dxiivmr + {ve{x)n (1-32) 

Finally, incorporating the interaction terms into the free Hamiltonian 
Eq. fOnl) one gets 

H=^ I dx[-{\j(t){x)f + UK{xje{x))% (1.33) 

where u is defined through u = vp^y{l + gi/2T[VFY — {g2/'2'^VFY and has 
units of velocity, while k is a dimensionless parameter which is defined by 



A+94/2.VF-92/2.VF_ 
^' 1+ gj2%VF + g2/2'KVF 



''Historically the interaction processes were given the notations 51,32, <73 and 54. How- 
ever, gi has the same form, for spinless electrons, as 52, and 53 denotes umklapp processes, 
which are not included in the Luttinger model. 
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The value of k describes the interactions in the system, k > 1 for systems 
with attractive interactions, while k < 1 for repulsive interactions. The value 
K = 1 corresponds to a non-interacting system. 

In conclusion, we see that a free-bosonic Hamiltonian can be used to de- 
scribe the physics of an interacting ID system. Thus the physical properties 
of such a system can be exactly calculated. The next significant milestone 
was by Haldane in 1981 [12], when he has proven that Luttinger's theory 
defines a new universality class which can be used also to solve other sys- 
tems. The problems which belong to this universality class, and for which 
the Tomonaga-Luttinger theory can be used, are called, after Haldane's sug- 
gestion, Luttinger liquids (or Tomonaga-Luttinger liquids). 

1.4 Motivation and Outline 

In the previous sections we have briefly described a few physical properties of 
ID lattices and of QDs. In the following chapters of this thesis we will delve 
deeper into these two systems, demonstrating some new findings for both, as 
well as discuss new features which arise once we couple both systems to each 
other. 

Recently, in the work of Kane and Fisher [ISj, it was shown that by 
coupling an impurity between two Luttinger-hquid leads the conductivity 
vanishes. Namely, the Luttinger liquid has a metallic behavior only in its pure 
form, and once disorder is introduced it immediately becomes an insulator. 
Indeed, it was recently demonstrated that the conductance through a ID 
disordered device produces infinitely sharp Coulomb blockade peaks at zero 
temperature [41] . 

So can we say that nothing happens when a QD is coupled to a lead which 
is described by the Luttinger liquid theory? Clearly, as we just said, the con- 
ductivity should vanish. Yet, one can ask what happens to other properties, 
such as the broadening of the dot's level, and the charge distribution in the 
lead. 

Does the coupling leave a signature on the dot's energy levels? Can we 
expect to measure a broadening of the levels as in a regular weak-coupling 
regime, or maybe the discrete form of the levels should re-appear because 
of the Kane-Fisher result? What happens to the charge distribution inside 
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the lead when we connect the QD? Does it show Friedel Oscillations? What 
happens when the lead is disordered? What happens to these phenomena 
when the lead's phase is changed between Tomonaga-Luttinger liquid and, 
let's say, charge density wave? 

In this work we try to answer these questions. In chapters [3] and H] we 
will examine a ID lead of spinless fermions which is connected to a QD. We 
restrict the QD, for the sake of simplicity, to have only a single energy level. 
While in chapter [3] we deal with a clean system, in chapter H] we investigate 
the influence of a weak disorder on such a system. In chapters [5] and [6] we 
ignore the leads, and examine closely two kinds of properties related to an 
isolated QD. In chapter [5] we present the results of an innovative numerical 
method for the ground-state energy of a disordered QD. We also utilize this 
method in order to examine the dot's addition spectrum, comparing our 
results to those obtained by the Hartree-Fock approximation. Chapter [6] is 
dedicated to the investigation of the g-factor of doubly-occupied QDs, in the 
presence of spin-orbit coupling and electron-electron interactions. 

Since much of the work utilizes numerical calculations, we choose to begin, 
in chapter [21 with a detailed review of the arsenal of numerical methods which 
are to be used. 



Chapter 2 
Numerical Methods 



2.1 Introduction 

In this chapter we present the methods with which the numerical resuhs of 
the following chapters are obtained. The main physical quantity which we are 
generally looking for, is the energy spectrum of our systems, including wave 
functions and their energies. Four methods are detailed here, and while their 
basic goal is quite the same, they have very different assumptions, algorithms 
and frameworks, and thus, they are utilized in different situations. 

In a few simulations the lowest energy state - the ground state - will be 
sufficient, while in other cases we might need to calculate the entire spectrum. 
In some of the systems we will be able to find exactly the physical properties, 
while in others we are forced to use approximations. And finally, some of the 
algorithms use the real-space as their framework, while others are built upon 
the momentum-space. 

In the following sections we describe in detail the four methods utilized. 
We start with the good old "exact diagonalization" method, which serves 
very successfully for the last few decades. Following it we present another 
very well-known method, the Hartree-Fock method, with which interactions 
between the electrons can be approximately treated. We then present a newer 
method, the density-matrix renormalization-group method, which was first 
established by S. R. White in 1992 [IS, |16], and which since then has got 
several hundreds of applications by a few dozens ofgroups around the world. 
We finish by presenting a bit more sophisticateqj version of this method. 



^and much more complicated 
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the particle-hole density-matrix renormalizat ion-group, which was used only 
rarely, but seems a very promising method. 

2.2 Exact Diagonalization 

The most simple method in which a system's energy spectrum can be cal- 
culated is by an exact diagonalization of the system's Hamiltonian. Numer- 
ically, the Hamiltonian is represented by a matrix, and the physical Hamil- 
tonian diagonalization is nothing else but a mathematical matrix diagonal- 
ization. In order to write the physical Hamiltonian, and to represent it by a 
matrix, one can choose the basis in which physical states are written, and ac- 
cordingly the way in which the matrix indices will be related to these physical 
states. 

Let's start with a simple example, that of a tight-binding Hamiltonian of a 
two-dimensional (2D) lattice with A rows and B columns, which is commonly 
used in order to represent a lateral quantum dot (QD) (see chapter [1]). We 
denote the total number of sites by (i.e., = AB). The motion of the 
electrons is by hopping from one lattice site to one of its nearest neighbors 
(NNs), which we denote as usual by (■■■). For simplicity we assume here 
that the QD is clean (without disorder), and we set the zero energy level 
at the constant energy of the electrons when they are located at the lattice 
sites. An intuitive choice of basis states leads to a notation in which state 
l^) refers to an electron in the i-th site of the lattice, so that 1 < i < A^. 
Denoting by Oj and aj the physical operators for creation and annihilation 
of an electron at site j, one gets the Hamiltonian of the system (compare to 



where t is the overlap integral between NN sites, and which is assumed to be 
identical for all sites. 

We now move on to writing down the matrix related to this Hamiltonian. 
Using the same indices notation, we will use the vector indices 1, 2, A^ to 
represent the states |2), |A^) and the matrix element {i,j) will now 
represent the physical quantity {i\H\j). The size of the Hamiltonian matrix 



Eq. dnD) as 




(2.1) 
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is thus M — N, and diagonalizing the Hamiltonian is done by a diagonal- 
ization of the M x M Hamiltonian matrix. The results, the eigenvalues and 
eigenvectors we get, are the energies and the wave functions of the system. 

How far can one go with such a method? A computational limit may 
result from the fact that an M x M matrix diagonalization is an O(M^) 
process, so that if, for example, the matrix size is doubled, e.g. for treatment 
of larger systems, or in order to consider the spin of the electron, the com- 
putational time is multiplied by a factor of 8. Nevertheless, with the current 
technology, the diagonalization task is still possible for quite large physical 
systems. 

Problems begin when interactions between electrons are considered. The 
process we have detailed so far, can help one find only the single-particle 
energies and wave functions. Let's assume that we'd like to put 2 electrons 
in a lattice of A rows and B columns. If there is no interaction between the 
electrons, then the many-body wave function is a Slater determinant of the 
single wave functions. In other words, when we put the first electron in a 
single-level \i) and the second in the resulted wave function is — 
73 — with an energy Eij — Ei+Ej. But when the electrons have 

some kind of an interaction between them, this is completely wrong. In that 
case the many-body wave function is more complicated, and in order to get 
the correct wave functions for two interacting electrons one has to include all 
the states in the matrix indices. In general, for Ue interacting electrons 
with spin, in a lattice of A rows and B columns, the matrix size will be 
M — (^^^) , with an exponential growth as a function of the system size, and 
it becomes too much for our simple 0{M^) diagonalization method, even for 
modest system sizes. 

2.2.1 Lanczos Diagonalization Method 

A partial solution for that problem is achieved by changing the diagonaliza- 
tion algorithm we use. If we agree to pay some price, we can use a more 
efficient diagonalization method, which works in about O(M^) steps. The 
price we pay is that we can't get the entire spectrum with such a method, 
but only a few of the levels. The profit is a decrease in the calculation time, 
which results in an increase of the size limitation for our systems. 
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There are two such diagonahzation algorithms that are more efficient 
than the 0{M^) one, the Lanczos method and the Davidson method. Both 
methods retrieve only a few eigenvectors, the ones with either the highest or 
the lowest eigenvaluej^. Both also assume that the matrices are sparse, and 
require an external efficient matrix-vector multiplication procedure. 

In the numerical work reported in this thesis, the Lanczos algorithm is 
used |47j . The idea is to replace the original matrix that we want to diagonal- 
ize by a new tridiagonal matrix with the same eigenvalues and eigenvectors, 
so that the diagonahzation task will be much easier. Furthermore, since one 
usually needs only the lowest or highest eigenstates, the size of the tridi- 
agonal matrix T can be much smaller than the original matrix size. One 
uses a recursive relation which is known as " Lanczos recursion" , in order to 
produce the new set of vectors ("Lanczos vectors"), with whom the original 
matrix transforms to the tridiagonal form. The Lanczos vectors are orthog- 
onalized using a Gram-Schmidt process, and during their production, a few 
tridiagonal matrices T, of varying sizes, are kept. These T matrices are then 
diagonalized iteratively, starting from the smallest, and the eigenvalues and 
eigenvectors obtained are mapped onto those of the original matrix, which is 
usually larger. If the mapping is not successful, a new iteration starts with 
the next T matrix that was kept [H]. 

To conclude, by using the Lanczos method one can enlarge the system 
size and try to include also interaction between electrons. However, there is 
still a serious size limitation when using exact diagonahzation methods, and 
the largest interacting systems investigated so far with that method were of 
6x6 with 4 spinless electrons. Of course that in special Hamiltonians with 
unique symmetry circumstances one might be able to increase this limit, yet, 
in regular cases the value of the upper size limit is very disappointing. 

2.3 The Self-Consistent Hartree-Fock Method 

The Hartree-Fock (HF) approximation is a well known method to deal with 
interactions in a "mean-field" way. The idea is to replace the exact potential 
felt by an electron due to all the others with a more convenient term. Instead 

^Actually they can retrieve the entire spectrum, but then they won't be efficient com- 
paring to a regular exact diagonahzation. 
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of solving the exact Hamiltonian and find the real many-body wave functions, 
one uses the HF approximation in order to find a basis of "one-body" wave 
functions which incorporate implicitly the interactions. 

Suppose for example that we want to calculate the spectrum of a QD, 
which is represented as an interacting 2D system of spinless electrons without 
disorder. Such a system is governed by the Hamiltonian (see chapter [T]) 

H = -t"^ {dl^an + H.c.) + ^ Va)^diand^, (2.2) 

(m,n) (m,7i) 

where the first term, to be denoted by Hq, represents the hopping matrix ele- 
ments between sites m and n which are NNs, and the second {Hint) represents 
NN interactionjfl. 

Writing the interaction term in k-space by using the relation dm = Ylia ^a'Paijn)) 
where is the a-th wave function (a runs over the momentum-space levels), 
gives 

^ = ^0 + ^ E 0:m0;(^)07H05M- (2.3) 

Instead of rewriting the entire Hamiltonian with {flf^ indices, Hartree has 
proposed that when we are interested on site m, we'll do an average, in the 
vacuum state, over the indices of the other site, i.e., take (0| h^p(j)*p{n)h^(j)^{n) |0). 
The meaning of the vacuum state |0) is a summation over all the states below 
kp, so that (0| h\(p*p{n)h^(t)^{n) |0) = Y.k<kp'Pl{'^)'Pk{n). We thus get 

^^r^'^ = EE^i'^^M^^'^^M^E'^^^H'^'^H (2-4) 

{m,n) a,5 k<kp 

= dl^dm V Y <Pkin)Mn), 

{in,n) k<kp 

which tells that one should add V^^^^n) "l^kKkp ^*ki.^)4>kin) to the matrix 
element H {171,171). 

A similar argument, which was proposed by Fock, leads to an addition of 
the non-diagonal elements, but with a minus sign because of the commutation 
relations. The averaging is now done over (0| b^f^(f)*p{n)bs(t>g{m) |0), thus giving 

^The modification of tlie interactions to Coulombic form is straightforward. 
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Efc<fc^ ^l{n)Mm), so that 

HfnT" = $^$^6l.C(m)6,0,H\/5^</)^H0,(m) (2.5) 



and therefore we should subtract V^f^^f^^ (l)l{n)(j)k{m) from the matrix ele- 
ment H{m,n), where m and n are nearest neighbors. 

The Hamiltonian matrix is now written in the A ■ B real-space indices. 
We first copy all the elements of Hq which are easily calculated in real-space 
coordinates. For the interaction we have these two contributions of Hartree 
and Fock terms. In the diagonal elements we should add the Hartree term, 
given by Eq. (12. 4p . and from the non-diagonal elements we should subtract 
the contribution of the Fock term, Eq. (12.51) . But unfortunately the 0's, 
which are required for these two terms, are the unknown wave functions we 
seek for. 

If we could have written this Hamiltonian, then its diagonalization would 
give us the correct wave functions. So what we do is to employ a self- 
consistent method: we take some initial wave functions (it can be simply the 
one-body wave functions), substitute them in Eqs. (12.41) and (12. 5p . and then 
diagonalize H and find a new set of wave functions. We then take this new 
set and calculate again Hint and H, and diagonalize H again. These steps 
are repeated until we get a stable solution, i.e., the output wave functions 
are identical (up to a defined accuracy) to the input wave functions. 

At the end, after the convergence of the self-consistent method, we have 
a new set of "one-body" wave functions which consider the interactions in a 
mean-field way. With this new basis one can compute the physical properties 
in a simple way. For example, since this is effectively a "one-body" basis, 
the ground-state energy of He electrons occupying the lattice is related to 
the sum of the lowest eigenenergies in the new basis. However, it is not 
identical to this sum, since one should compensate for double counting in the 
Hartree and Fock terms. It is easy to prove that the ground-state energy is 



{m.n) a, 7 



k<kp 
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thus 

Egs = (2.6) 

k<kp 

{m,n) l,k<kp 

where the first sum runs over the single-particle eigenvalues, and the second 
one counts for the Hartree and the Fock terms [lOj . 



2.4 The Density-Matrix Renormalization Group 
Method 

In this section we describe a method which is used in order to get the low- 
est lying levels of a system without the usage of exact diagonalization, the 
density-matrix renormalizat ion-group (DMRG) method. This method works 
for much larger systems than those one can treat by using the exact methods, 
and gives much more accurate results than those one can get by using the 
HF method. 

The DMRG method gives only an approximate result for the physical 
properties. However, it has proven itself as very accurate in various systems. 
Nevertheless, when one starts an investigation of a new system, they must 
assure that the results are accurate and represent the real physics. Another 
severe limitation of the DMRG method is its implicit restriction to one- 
dimensional (ID) and quasi ID systems. 

The DMRG method is actually a development of some earlier numerical 
renormalization group methods. The basic idea behind many of the nu- 
merical renormalization group (NRG) methods is very simple: We cannot 
diagonalize the full Hamiltonian since it is too large. However, if we need 
only the lowest levels, and if we can increase the Hilbert space iteratively (for 
example by adding a site after another) we can omit in each step the states 
which seem the most unimportant, and then continue with just part of the 
Hilbert space. Afterwards we add another site, the Hilbert space increases 
again, and we again throw the unimportant states. We continue this process 
site after a site, until we finally get to the size of system we want. If we 
did a good job during the iteration process, i.e., all the states we've omitted 
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had indeed a negligible contribution to the final ground state, then we have 
a good approximation of the ground state. 

So the million-dollar question is how to decide which of the states we'd 
like to keep, and which can be omitted. The old methods used a pretty logical 
criteria, in which the states being kept are those with the lowest energies. In 
fact, this method was first developed by Wilson for the problem of an 
impurity coupled to a ID non-interacting lead. It was shown that the system 
can be mapped into a new system with a tridiagonal form, representing an 
impurity level which is coupled either to nearby sites, or to states with very 
low energies. In addition, it was proven that the hopping elements of the new 
system decrease algebraically between consecutive sites. In such a problem, it 
was shown that the NRG method gives accurate results. However, attempts 
to use this method for cases in which the leads contain non-trivial physics 
(for example when interactions between the electrons in the lead are not 
neglected) weren't so fruitful. 

One of the reasons for the NRG failure in those cases, is that there is no 
guarantee that a state which wasn't important for a 10-site system, is not 
important also for the 100-site one. Furthermore, the boundary conditions 
are changed during the process of sites addition. The states of the A^-site- 
system iteration do not have the same boundary conditions as those of the 
iteration with N + 1 sites. 

2.4.1 DMRG: Quick Overview 

In 1992, S. R. White proposed a different approach regarding the decision 
which states to keep [ISIHB], the DMRG method. Since then, it has become 
one of the most useful numerical methods for treatment of ID systems, with 
many applications and variations [50| |5T]. 

The main idea of this method is to include in the Hamiltonian of the 
current iteration, not only the "system" sites that were already added, but 
also some "environment" sites which are coupled to the current system. Not 
only does the inclusion of more sites assure the correct boundary conditions 
at each iteration, but it also forces the diagonalization process to take into 
account information from larger system sizes. This leads to the diagonal- 
ization of a Hamiltonian representing a larger physical system, thus getting 
a ground state of a larger basis. To return to the correct system size, one 
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The Superblock 



System Block 


Environment Block 















Figure 2.1: The DMRG superblock composition. The system block (in 
black) , containing the dot (the most left site) coupled to a lead with — 1 
sites, is coupled to a new site (in blue). Then the entire system block is 
coupled (green line) to its mirror picture (red sites) which represents the 
environment to form the superblock. For a block of p states, one gets a 
superblock of states. 



now defines a density matrix from which tlie environment degrees of freedom 
are traced out. By a diagonalization of this density matrix one can identify 
which of the system states are the most important for the following steps, 
since they get the highest eigenvalues. 

Before getting deeper into the technical details of the method, let's de- 
scribe in general the iteration process. We start an iteration when our system 
already contains a block of sites, whose physical operators are represented 
by matrices of size p. This is named the " System Block" schematically shown 
(in black) in Fig. 12.11 We assume that this size, p, is the constant number of 
states we want to keep during the iteration procesqj. Now we would like to 
add the next site (the blue one). We thus enlarge our basis set by a factor 
of 2 (assuming for the moment that the electron is spinless), resulting in an 
enlargement of all the operator matrices to the size 2p. The idea of White 
is, as mentioned above, to diagonalize a larger system's Hamiltonian, called 
a " superblock" , by an addition of some environment sites. This can be done 
simply by taking a mirror picture of the system we have as the environment 
(the red sites), connecting the system and environment by the regular Hamil- 
tonian terms (represented by the green line). This results, in principle, in a 
superblock size of 4p^, which is now diagonalized (usually using Lanczos or 
Davidson methods), in order to find its 4p^-sized ground state. From this 
ground state one constructs a density matrix of size 2p, by summing over 



^Increasing the value of p will increase the DMRG accuracy on one hand, but on the 
other hand it affects dramatically the required calculation time and memory resources. 
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the environment indices, so that the remaining 2p indices represent our ini- 
tial block with the addition of the new site. By taking the upper half of 
the density-matrix eigenvectors (i.e., those with the larger eigenvalues) as a 
projection operator, one can now transform all the 2p-sized operators to the 
new basis, of size p. This new block is the starting point of a new iteration. 

2.4.2 DMRG: a Technical Overview 

We now move to describe the technical details related to the Hamiltonian that 
we use in different parts of the thesis. This Hamiltonian describes the motion 
of spinless electrons on a ID lattice ("lead") of size L, which is connected at 
one side to a QD. Such a Hamiltonian should consist of both Hqd and Hiead 
described in the previous chapter (Eqs. (11.11) and (I1.13P ). We now restrict 
the QD to have only a single site, and in addition, we add a hopping term 
between the dot and the first site in the lead, whose matrix element will be 
denoted by V'o, If no interaction . assumecg one gets 



J = l 

where and a represent the creation and annihilation operators of an elec- 
tron in the single state of the dot, of energy eo. We now use the DMRG 
method in order to get the ground state of the Hamiltonian Eq. (12.71) . 

Step 1: Building the First Block 

In the very first step one treats a very short system, e.g. a dot level coupled 
to a lead with only one site. Since this system has in total 2 sites, it has only 
4 basis states. The 4x4 Hamiltonian matrix of this small system is now 
written exactly: 



We now want add a new lead site to the system, thus increasing the 

matrix size by a factor of 2. Actually we can assume that we are in the 

beginning of a new general iteration, since the next steps do not depend on 

the number of sites the system already contains, specifically this is of course 

true also for the first iteration. 

^The inclusion of NN interaction is straightforward, see the discussion in the next 
subsection. 



L-l 




(2.7) 
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Step 2: Addition of a New Site 

Let us denote the current block matrix size, containing already n lattice sites, 
by p. This block is represented by the black sites in Fig. 12.11 We now add 
the next site, number n + 1 (the blue site), and the new matrix size after the 
addition is 2p. The new Hamiltonian, denoted as Hb,, should now contain 
the same elements as Hb, the diagonal elements of the new site, and also the 
connection between the new site to the previous one: 

Hb, = Hb- i(4cn+i + cl+iC„). (2.9) 

Since the expression of the Hamiltonian is done using the operators cj^ 
and Cn, one can see that there is a need to store the operators of the previous 
step. Later we shall see that the operator cj^c^ may be also used during the 
calculation. Since the DMRG process change the basis states continuously, 
and the states kept are not a complete basis of the Hilbert space, in general 

) ^ ^{i\ci\m) {m\cn\j) , (2.10) 

m 

therefore this operator, as well as every other operator that will be required 
later, should be specifically kept. 

Aside from increasing the size of H, we perform a similar procedure for 
any operator we shall need. These operators are rewritten in the new 2p-sized 
basis, which includes the previous basis and the new added site. In the next 
steps we will, in general, decrease the Hilbert-space size back to the original 
size p% 

Step 3: The Superblock 

We now have a system of n+1 sites. Recalling that the DMRG idea is to place 
this system inside some sort of an environment, one can see that the simplest 
way to do it is to take a mirrorpicture of the system as the environment, 
and couple these two subsystemq^. We thus refiect the B» system in order to 

^At the first few iterations of the process, the size of the new basis might be still smaller 
than p, the matrix size we want to keep. In this case we continue to add sites, each time 
doubling our matrices size, until the matrix size is greater than p. Then we continue to 
the next two steps, which are needed for the truncation of the Hilbert space. 

^Alternatively one can grow a clean lead without a dot at an initial DMRG stage, and 
then couple it to the Bu block at its right edge. 
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get an environment block ^B^ {B^ denotes the reflected block), and couple 
the two new sites in the middle by the regular hopping term (see Fig 12. II) . 
By that we get the superblock Hamiltonian 

of size {2pY = Ap^. 

Step 4: The Density Matrix and Hilbert-Space Truncation 

The next step is to find the ground state \E'o of the superblock, usually by 
using the Lanczos algorithm!^. Dividing the indices of the ground state to 
\E'o = Tl,is jE^^sjE\'^s)\3E) 1 where S (E) denoted a system (environment) 
index, we use to build the density matrix p by tracing out the environment 
degrees of freedom: 

3E 

The density matrix, whose size is 2p, is now diagonalized and the p eigen- 
vectors which have the highest eigenvalues are used as a projection operator 
O of size 2p X p. If we denote p\ua) = Ua\ua), then we choose the eigenvec- 
tors \ua) with the highest UaS to form an 2p x p matrix O. For every stored 
operator Q we now use O^QO ^ Q in order to decrease its size to the initial 
size, as it was at the beginning of the iteration. 

Now that the system block has n + 1 sites, and all of the operators are 
back in their initial size, we can return to step (2), and add the next site. 

As steps (2) — (4) are repeated iteratively, the physical size of the system 
increases, while the Hilbert-space size remains the same. The iteration pro- 
cess is stopped when the physical observable we want to measure, e.g. the 
occupation of the dot level, converges. 

^Usually the superblock is not diagonalized as a 4p^-size matrix, but in sectors con- 
taining different numbers of electrons. This can be done since the Hamiltonian does not 
couple states with different number of electrons. Using such a method, one should diag- 
onalize each sector, and finally choose the eigenvector with the lowest eigenvalue among 
all sectors. 



2.4. THE DENSITY-MATRIX RENORMALIZATION GRO UP METHOD 35 



2.4.3 More Complicated Models 

The generalization of the previous subsection to the case of interacting elec- 
trons depends crucially on the details of the interaction term, specifically on 
its range. The inclusion of NN interactions is straightforward, since in each 
iteration we just have to add the term /cJ^Cncj^+iCn+i- The operators cl and 

are of the previous iteration, and since they are needed for the hopping 
term, they are kept anyway. The operators c]^^^ and c„+i are of the new 
added site, so they can be written explicitly. 

The same argument holds for the superblock creation, in which interaction 
should just be added between the two new sites (in the green line of Fig. 12.11) . 




Figure 2.2: The DMRG superblock composition for next-nearest- 
neighbor interaction model. The terms added in step 2 are drawn in orange 
(t, I), blue and violet {I2), while the terms of step 3 are colored green (t, /), 
magenta and maroon (12) • 



If the interactions include terms of longer range, the case is different. 
Let's look at the model of next-nearest-neighbor interactions, in which there 
are two interaction terms, one for sites which are NNs (denoted by /) and 
the second {I2) for sites in distance of 2 lattice sites. 

The simplest solution to that case is that in every iteration we add two 
sites, instead of a single site as explained in the previous section. In addition, 
we will store in memory the operators of the two sites added in the previous 
iteration. Let's denote the new added sites by Ni and N2, and the sites 
added in the previous iteration (and whose operators are kept) by Pi and P2. 
A schematic picture of the superblock is shown in Fig. 12. 2[ In step (2) the 
addition of the new sites will be simple, where we will make the following 
connections: A^^i will be connected to Pi with I2 and to P2 with J, N2 will be 
connected to P2 with I2, and A'^i and N2 are also connected with I. 

In step (3) the superblock formation is also as simple as that, since there 
are 4 new states in the middle, and the interactions do not connect any other 
sites. So actually the algorithm is identical to the one detailed above. The 
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main difference is the sizes: the addition of 2 new sites increases the operators 
basis by a factor of 4, so that the superblock Hamiltonian size, which is to 
be diagonahzed, is IGp^. This decreases, in principle, the upper limit of p, so 
that the accuracy might decrease. 

Such a solution works, of course, only for a limited range of interactions. 
If we want to include more interaction terms, it quickly becomes impossible 
to increase the number of sites we add in each iteration. Another solution 
might be to add in each iteration only a single site, but to keep the operators 
of all previous added sites in memory, and use them in the connection of the 
new site, and in the formation of the superblock. However, in this solution 
steps (2) and (3) require much heavier calculations, e.g. multiplications of 
large matrices, and they also include mixing of operators from the system and 
the environment sides. In general, this leads to a much slower calculation. 

As a last remark we note that the simple solution we described for the 
next-nearest-neighbor interactions model, i.e., by adding 2 sites at a time, is 
practically similar to a solution for the case of spin 1/2 electrons with NN 
interactions, in which instead of adding two sites, we add in each iteration a 
new single site with 2 spin states. 

2.4.4 The Finite-Size DMRG 

While the accuracy of the DMRG method described above was pretty good 
for several systems, there are other problems for which the DMRG results are 
not good enough. In 1993, only a short time after his first publication of the 
DMRG method. White has proposed an improved version of the algorithm, 
the finite-size DMRG [56]. 

The idea is to use the DMRG process on a fixed-sized system, instead of 
an infinite one. The application is based on the DMRG iteration process, 
but requires more computer resources, mainly memory. In this method the 
operators of every iteration, i.e., of each and every size of the physical system 
(e.g. the lead), are stored for future use. Let's assume we'd like a system of 
a single- level dot coupled to a lead of size L — 1, so that the total length of 
the system is simply L. During the iterations, one can combine two saved 
operators by connecting operators of the left hand side of the system of size 
X to those of the right hand side of the system, with size L — x. Enumeration 
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over X from x = ltox = L — 1 will give us a complete path through the 
lead. 



(a) ®# O 



(b) ® 



(c) ®0 



o 



(d) ® 



Figure 2.3: A sweep of the finite-size DMRG algorithm, (a) Starting 
from the minimal system size, (b) the system grows and the environment 
size decreases, (c) the minimal size of the environment was reached, (d) the 
environment growth. The entire process of (a)-(b)-(c)-(d)-(a) is a sweep. 



Let's now explain exactly how it works. At first, we build a lead of size 
L by the regular DMRG method. The only difference is that in the process 
we keep a copy of all the operators representing the physical systems of sizes 
1, 2, L. We denote these operators as the environment operators, and they 
represent the right hand side of our desired system - see Fig. 12.31 - a clean lead. 
The operators which represent the left hand side of the system, containing 
the dot and the left part of the lead are called the system operators. 

We now start again with a new system from the left hand side, the part 
of the dot, by taking a new first block composed of a dot and a single lead 
site. This is the first block, whose physical length is 2. To this block we add 
a new lead site, and now we'd like to get a superblock. In the regular method 
we would take this 3-sites system and couple it to its mirror picture, getting 
a superblock of length 6. Instead, we couple it to a block which will make 
the total size exactly L. Such a block is composed of a previously-stored 
environment-block of size L — 4, and an extra new site. The total physical 
length of the system which is diagonalized - the superblock - is thus L. A 
schematic picture of this process is given in Fig. 12.3( a). 

The superblock is now diagonalized and the density matrix is formed. 
The projection process over half the space gives the system operators for a 
3-sites system. We now continue by adding site #4, and coupling the new 
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system to an environment of size L — 5, in a similar way, and finally getting 
the system operators for the next iteration. We can continue this process 
iteratively (see Fig. 12.3( b)) until we reach the other end - when our system 
has L — 4 sites on its own, and the size of the environment is the minimal 
size - 2 sites (Fig. EJl^c)). 

During the process we keep copies of all the system operators, since we are 
just about to use them: after reaching the right end, we want to come back 
from right to left (see Fig. 12.3( d)). connecting the new growing environment 
to the stored system operators, while the system side goes back from size 
L — 4 to 2. During these steps we update the environment operators, to be 
used in the next time, when we'll move from left to right again. 

We can continue this process in a zigzag manner, going from the dot side 
to the right, and then back all the way, where each back and forth iteration 
- (a)-(b)-(c)-(d)-(a) in Fig. 12.31 - is called a sweep. When we are moving to 
the right we update the system operators, and when we are on the way back 
we update the environment ones. 

Usually after a few sweeps the physical observables converge, the results 
are sufficiently accurate, and then we stop the iteration. 

2.4.5 Finite or Infinite 

The last subsection about the DMRG method will be devoted to the number 
of particles (ue) in the "finite-sized" system. Most DMRG applications use 
TIq clS db quantum number which is fixed during the renormalization process. 
In practice, one can order the block states in a block-diagonal form by their 
electrons number, and conserve this form when the superblock Hamiltonian 
is created. Furthermore, one might restrict the superblock state to have a 
certain number of electrons, by coupling only states which give together the 
required number. In such a way, and since the Hamiltonian does not couple 
states with a different number of particles, will remain a good quantum 
number even after the truncation. This of course has some impact on the 
superblock composition, which now has a block-diagonal form where each 
block contains states with a certain number of particles. In general, it makes 
the numerical method work much better, though it is mainly relevant to cases 
in which there is an external constraint on ng. 
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In our DMRG application, however, we intentionally do not keep He as 
a constant. Instead, we diagonalize the term H — fiNe, where Ne is the 
number operator, and fx is the chemical potential. We keep constant 
and let the renormalization physics determine what ng will be. In practice, 
we form superblock states of different particle numbers, and we do not keep 
the order of states in the blocks, and this has an important impact on the 
renormalization. Even though in the superblock ground state of the first 
iteration we get an integral number of particles, it does not need to remain 
integral after the truncation process takes place, so that in the next iterations 
the number is not a good quantum number any more. 

A direct advantage of this method is the option to get non-integral values 
for rie- As will be shown in the following chapters of this thesis, this prop- 
erty can be essential for obtaining, from a finite system, results which are 
accurate for an infinite system. As a simple example, let's think of a phys- 
ical system which has a particle-hole symmetry, so that with taking /i = 
should be exactly half filled. Without interactions, such a system can be 
exactly diagonalized (even if it is a large one), and the half filling property 
can be exactly checked. However, if the size of the system L is odd, even for 
L ^ 1, an exact diagonalization cannot give the "correct" answer (which is 
obviously correct just for the infinite case). With our implementation such 
an "infinite-system result" can be reproduced even by using systems of only 
a few hundred sites. From the experimental point of view, the results we 
obtain with this method, for example on ID systems, can effectively describe 
a situation of a finite section of a ID wire which is coupled to an external 
electron reservoir. 

2.5 The Particle-Hole DMRG Method 

The main disadvantage of the DMRG method is the need to have a well- 
defined order in which the system is enlarged. Regularly this limits the use 
of DMRG to ID, or at least quasi ID, systems. 

In fact, a short time after the first appearance of the DMRG, a different 
scheme of the DMRG process usage was proposed [52]. If one adopts the 
momentum-space point of view, instead of that of the real-space, a new 
meaning replaces the traditional ID approach: the order of the energy states. 
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The order in which states are added can now be defined physically by their 
energies, leading to the opportunity to start from the lowest lying states, and 
add the next ones iteratively. In each iteration one adds a new energy state, 
thus increasing the size of the Hilbert space, and then decreases this size to 
the original size, in a similar fashion to that is done in the regular DMRG 
process. 

Since the first idea to use the momentum-space was proposed in Ref. [52] . 
some implementations of this method were written in various fields [531 \5^ . 
In particular interest for us is the implementation named particle-hole DMRG 
(PH-DMRG), since we utilize in this thesis a similar idea for QDs. Sev- 
eral implementations of this sophisticated method have been reported so far, 
focusing on ultrasmall superconducting grains [551 IMl EH |58] and nuclear 
shell- model calculations [59l [601 El] • Some of these results, accompanied by 
a brief introduction to the algorithm, can be found in [621 [63] . 

2.5.1 PH-DMRG: Quick Overview 

The name particle-hole DMRG is based on a special implementation of the 
momentum space, where the renormalization starts from the Fermi energy 
and sweeps to both sides. Lets assume that we want to find the ground state 
of a system composed of a 2D lattice of A rows and B columns, in which 
there are Ue interacting spinless electronqj. An exact diagonalization of such 
a system requires, as noted above, the diagonalization of an [^^) matrix, 
which is often too large. 

So now let us describe how the PH-DMRG proceeds. First, a single 
electron basis should be chosen, and the most intuitive choice is the sin- 
gle electron eigenvalues, arranged according to their energies. However, the 
single-particle eigenvectors do not consider the interaction at all, and it can 
have a huge impact. A better way is to use the Hartree-Fock eigenvectors as 
the starting point to the PH-DMRG process. 

So now we have a set of wave functions sorted by their energies. In 
order to take into account exactly rig particles, the simplest way is to fill the 
lowest He states with electrons, leaving all the other states empty. This is 
depicted schematically in Fig. 12.4( a). The idea of the PH-DMRG process 
is to change this basis of states, in order to decrease the energy. This is 



^The choice of this example here is based on the usage of this method in this thesis. 
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done by considering some higher states also, while controlling the number of 
particles. Namely, putting a particle in a new energy state requires taking 
out an electron from a filled level, or, in other words, introducing a hole in 
one of the low states. 



P States 
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(b) 



(c) 



Figure 2.4: The PH-DMRG superblock composition for the first itera- 
tions, (a) The energy states after employing the Hartree-Fock algorithm, 
showing the occupied states (below Ep) and the unoccupied ones, dividing 
the states to P- and H-blocks. (b) The first superblock composition, com- 
posed of the 4 levels around Ep. (c) The superblock composition of the 
next iteration, after basis change of the 4 levels from the previous iteration 
(dashed lines) and the addition of 2 new levels. 



Therefore we keep two lists of energy states. The first one contains the 
energy states above the Fermi energy, which are initially unoccupied. These 
states will be called the "Particle" states (or simply "P-states"), since one can 
add a particle to the system by filling a state from this list. The second list 
is of the "Hole" states ("H-states"), which are the (initially) occupied energy 
states. In the situation we've just described, there will be total of = Ue 
hole states, and = AB — rif. particle states. We maintain the "Particle" 
and the "Hole" lists separately during the process, creating two kinds of 
blocks, which are combined once in a iteration to form the superblock. 
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Now let's explain in general the iteration process. At first we start with 
a small number of energy states in the vicinity of the Fermi energy, for 
example we can take one state from each list to be in the relevant block, 
i.e., the particle block (P-block) will contain the first state above Ep and the 
hole block (H-block) the state just below it. This is of course true only in 
the first iteration, and we now want to present a general iteration process, 
so let's assume that we start a new iteration in which the number of energy 
states which are already contained in the blocks (whether it is the P-block 
or the H-block) are n — 1. 

We increase the size of the two blocks by adding a particle state to the 
particle block and a hole state to the hole block, thus increasing the number 
of states to n in each block. From these two blocks we build the superblock 
Hamiltonian by taking an equal number of states, /, from the two blocks, 
where / represents the number of particle-hole pairs. For example, in the first 
iteration (see Fig. 12.4( b)) there are 4 single-particle states in the superblock, 
let's denote them by with i = 1,2, 3, 4. There will be a single superblock 
state in which / = 0, i.e., without any electron- hole pair, so that the lowest 
two states are filled and the higher are empty, and this state can be denoted 
by 1 1,2). There will be four superblock states with one p-h pair: these are 
the states |1,3), |1,4), |2,3) or |2,4). And finally it includes the state |3,4) 
in which two p-h pairs exist. So in the first superblock, when 4 single-level 
states are treated, the Hamiltonian is of size 6, being (2). 

The diagonalization of the superblock can lead us towards a new basis, 
by using similar considerations as in the regular DMRG process. We project 
the superblock ground state over the space of the particle states by summing 
over the hole indices, and create a P density matrix. We do the same for 
the Hole states, by summing over the particle indices and get an H density 
matrix as well, and we use these two matrices in order to get a new P basis 
and a new H basis, which replace the ones that we've started the iteration 
with. In each state of the new basis we keep a copy of the number operator 
of particles (or holes) involved in this state, in order to be able to couple it 
later (in the next iteration) in the same manner. 

Then we are in a good position to move on to the next energy states. We 
increase the number of states involved in the game by taking the next states 
- a P-state will join the P-block, and an H-state to the H-block (Fig. 12.4( c)). 
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As the Hamiltonian is built out of operators, we should keep a copy of all 
the operators we need. In each step the number of energy states increases, 
and so is the number of required operators. Since every operator is an p x p 
matrix, the memory requirements can be a serious limit of the PH-DMRG 
method. 

2.5.2 PH-DMRG: a Technical Overview 

We now present the technical details in relation with the Hamiltonian we use 
in this thesis The Hamiltonian describes the motion of interacting spinless 
electrons in a 2D disordered lattice (a QD). We will write the Hamiltonian 

of the system as 

H = Ho + Hint, (2.13) 

where Hq includes the disorder and the hopping elements, and Hint is the 
interaction. In real-space we can write (see Eq. (11. ip ) 

i^o = emflL^m - t ^ {alna-n + H.C.), (2.14) 

m {m,n) 

where the first term contains the on-site energies, which are diagonal in the 
lattice sites (m and n denote both rows and columns for simplicity), and 
the second presents the hopping elements between sites m and n, where the 
notation (...) denotes NNs. 

The interaction term for NN interaction is (Eq. (11. 5p ) 

Hint = ^ al^alanam. (2.15) 

(m,n) 

In order to work in k-space, one needs to write all these terms in their 
k-space equivalence. Using the relation = Ylia^a4>a{^)-, one geti^ 

^e„aj„a^ = ^hihp^tn^(t)l{rn)(t)p{m)]. (2.16) 

m a,f3 m 



lOThe PH-DMRG method is used in chapter 5. 

-'^-'^The modifieation of the interactions to Coulombic form is straightforward, and will 
be discussed in chapter [Sj 

Along this section the operators a and the indices m, n are used for real-space, whereas 
the operators b and the Greek indices are used for k-space. The indices i,j,k,l will be 
used for the matrix indices. 
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The hopping term is less simple, but without too much effort one gets 
J2 («ln«n + aiaj = 2 J2 E C(^)</'/9(^), (2-17) 

(m,n) a,j3 (m,n) 

SO that we can write 

^0 -J2WT.'rn<l>l(m)<l>p(m) -2tJ2 (2-18) 

a,/3 m (m,n) 

Since shortly we'll begin multiplying the matrices which represent the 
operators and 6^ are p x p matrices) in order to get the Hamiltonian 
matrix, it is worth noticing that the term in the parentheses does not depend 
on the matrix indices, so it is calculated only once for each a and (3. 

In the same way it is easy to write the interaction element in k-space, as 

Hint = E biblb^bs V J2 K{rn)ct>%n)ct>,{n)Urn), (2.19) 

a,l3,^,5 {m,n) 

where again, the second summation is done once for each set of a, /3, 7, 5. 

In order to simplify future calculations we define the anti-symmetric in- 
teraction term 

- ^a{mWp{n) [.p^{n)Um) - h{n)4>^{m)] (2.20) 

(m,n) 

- <P*ain)(p*f3{m) [0^(n)05(m) - (j)s{n)(f)^{m)] 

(m,n) 

We then get 

H^nt = ^ ^ blblb^bs. (2.21) 

a,/3,7,<5 

After having the Hamiltonian written in our favorite notation, we can 
start describing the iteration process. The algorithm main steps are roughly 
similar to those of the regular DMRG method. Nevertheless, the details are 
very different, so we now explain in detail the iteration steps. Recalling that 
the first iteration is done after an application of the Hartree-Fock method, we 
assume the existence of a well-defined ordered set of wave functions (sorted 
by their energies). 
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Step 1: Building the First P- and H-Blocks 

The first blocks contain only one state each. Actually we create here the 
matrices for all the operators we'll later need. The Hamiltonian of the su- 
perblock contains only terms with an even number of operators, i.e., terms 
like ¥b and b^b^bb, but since we prepare the P- and H-operators separately 
while the Hamiltonian is built out of combinations of them, we must keep all 
the possible combinations, including those with an odd number of operators. 
We thus keep the following combinations: b, b'''b, bb, b'''b'''b and b'''b'''bb. The 
operators b^b^ and b^bb are omitted since they can be calculated using 
Hermitian adjoint relations. Every b in the above terms, in addition to rep- 
resenting ap X p matrix, should have a subscript Greek index, denoting the 
energy state to which it is related. Nevertheless, in the case of the four-state- 
operators one does not need to store every operator on its own, since only 
their summation — ^ ^ ^ Y^tal is required for the superblock. 

Therefore, only this sum, D4, is stored. In the very first iteration, there is 
only one P-state (ccq) and only one H-state (/3o), so that the list of operators 
we create is doubled: there are 5 P-block operators and 5 H-block operators. 

In practice, after building explicitly b and b^ , we can define the other 
ones by a multiplication. Note that this simple method works only in the 
first time, because after the first truncation the set of wave functions are not 
a complete basis any more, exactly as was the situation in the regular DMRG 
algorithm. Thus for example 

{i\ k 62 \j) ^ Y^{i\b^\m){m\b2\2), (2.22) 

SO we must keep all the operators explicitly in the memory. 
Step 2: Addition of a New State to the Block 

We now move to the next part: adding a new state. Let's think about adding 
a state to one of the blocks assuming there are already n — 1 states inside, 
and we are adding the state. What wc need to do is to add to the list of 
operators all the combinations with 6a„ or b\n- increase in this step 

the size of every operator by a factor of 2 in order to include the new state, 
so that the operator matrices are now of size 2p. 
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For simplicity, let's denote the operators of the previous basis (all the 6's 
and 6^'s of the first n — 1 states) simply by b, while ba will now represent the 
operator of the new state (instead of 6a„)- We do the same for the indices, 
so that the new state vector will be denoted by \i,ia), where i encodes all 
the "old" indices, and = 0, 1 is the index of the new energy state. We 
now wish to calculate the operator matrix elements in the new basis, and in 
general we will look at the element {i,ia\0\j,ja), where O can be any one of 
our operators. Since \i,ia) = \i)\'ia), we can write it as (ia|(z| O \j)\ia)- 

Let's start from the one-state operator b which is the easy onqlf|. If the 
operator O is of the old basis, denoted in general as b, one has obviously 

{^amb\J)\J^) = {^\b\J)5,^,^, (2.23) 

since b does not operate in the subspace of ia and ja- Note that {i\b\j) is 
already known (and saved) from previous iterations. On the other hand, if 
O = ba "we need at first to replace the ba and j , and since these are fermionic 
operators it gives a factor of (—1)"-', where rij is the number of particles in 
state j. Thus we get 

{^a\{^\ k \3)\3a) = {-IT^ 5,A^a\k\3a) . (2.24) 

where {ia\ba\ia) operates only in the subspace of the new energy state, so it 
can be written explicitly . 

It's about time to get to the two-state operators, hi) and })^b. Since the 
process is similar for both of them, we will discuss only the bb case, and it is 
divided to four possibilities: bjja-, kb, bba or bb. 

The idea is the same as in the one-state operators, so to make a long 
story short, we'll write here only the final expressions. Whenever there are 
two operators (or more) of the same basis we put an extra subscript, in order 
to differentiate between them. However, the order of operators is always 
conserved. 

{^a\{^\khma) = {Ahkm^,^ (2.25) 
(lambb^m^) = {-lp{^\b\J){^Ma) 

{iamkb\j)\jo.) = -{-lpmj){la\ba\ja) 
ba,l ba,2 \j)\ja) = &a,2 \ja)- 

^^As noted above, one does not need to keep b'' , since it can be calculated using b. 
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Note that the minus sign (in the third hne) results from the fermionic 
commutation relation between ba and b. 

Using the same method we can write down the 8 possibilities of the 
three-state operator bjj)l,bx in the enlarged basis. However, using the anti- 
symmetric fermionic relation y^J)\bx = —b^b^^Jbx for li ^ v (with bjfijfix = 0), 
one can keep only those operators with ^ < v. Therefore, the only relevant 
possibilities are 

{^a\{^\b\b\h\J)\Jo) = {A b\ bl k m^j^ (2.26) 
{^a\{^\b\bik\J)\J^) = -{-l)-^^\b\k\J){^a\bi\Ja) 
{^a\{^\b\blb^\J)\J^) = i-lp{t\blbl\j){tMja) 

{ia\{i\ h a,2|j)|ja) = (^1 

Similarly, from the 16 four-state operators y^Jb^bxb^ only those with yU < 
V and A < K are necessary, which again leads to only four possibilities. 
However, since we store only the sum of these four-state operators multiplied 
by V^i/Ak, which we have already denoted as D^, we first enlarge the basis for 
the case in which all of the operators are from the old basis, using 

{t^\{t\b,\j)\3^) = {t\D,m^,^. (2.27) 

Next we add to 1)4 all the new cases, in which u and / 01 k represent the 
new energy level. Namely, we add to {ia\{i\ -D4 

= «) (-l)"M^l^l^2&3|j)(^a|&J.|ja)V^W,3 (2.28) 

(«: = «) (-l)"^(^|&I&5^3|j)(^«|&a|j«)^l,2,3,a 
") - (^1 ^1 ^2 \j){ia\ y^al ba,2 \ja)Vi 

,01,2,02 5 



[V = K 



and finally we store only D4. 

At last, we remind that every operator b and listed in Eqs. f l2.23p - 
(12.281) should have an index which runs over all the states that have already 
entered to the game in the previous iterations. One should also note that 
such a block enlargement is done separately for both the particle and the 
hole blocks. 
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Step 3: The Superblock 

After the enlargement of the blocks we can combine the P- and H-blocks 
to form the superblock. In the superblock it's the time we really build the 
Hamiltonian of the system. This is done by matching for every P-state which 
contains / particles, an H-state with the same number / of holes. Since the 
operators of each block are matrices of size 2p, and they operate on different 
spaces, the maximal size of the superblock Hamiltonian matrix will be 4p^. 
The matrix indices are a direct product of a P-index and an H-index. 

So if we denote the P-indices by i for the row and k for the column, and 
the H-indices by j and /, accordingly, we now want to calculate the matrix 
element {j\{i\H\k)\l) . The Hamiltonian is given by = Hq + Hint, where 
Ho and Hint are written in k-space and given by Eqs. f l2.18p and (12.211) . in 
which the operators are either the P-operators or the H- ones. 

For the calculation of Hq we need to enumerate over the two state indices, 
a and /?, and divide the cases to four possibilities, depending on whether a 
and f3 are states from P- or H-blocks. By denoting the P (H) operators by 
an appropriate subscript we get 

{jmbip^pim) = {^\bipkpm,i (2.29) 
{mbUh,,\k)\i) = i-iT'^iti bip\k){j\bp,,\i) 
{jmbUh,p\k)\i) = -{-lrH^\h,p\k){J\bl,\l) 

{jmbl,kh\k)\l) = S.Aj\bUbf,,,\l). 

The calculation of Hint is just a bit more lengthy, because here we have 
4 state operators (a,/?, 7,5), which result, in principle, in 16 possibilities, 
depending on the configuration. However, because of the enumeration over 
these indices, and because of our definition of the anti-symmetric interaction 
term Vaf^-yS, some of the possibilities are equivalent. For example, let's look at 
the terms ppph (three P-operators and then an H-operator) and pphp. Using 
the fermionic relations and the antisymmetry of V, we can write 

K/376^L,p bl^p Kh hp = (-'^'a/357)(-^I,p b^ ^ bs,p Kh) , (2-30) 

and thus these two terms give the same contribution. In the same manner 
we get phpp ~ hppp (~ denotes an equivalence), as well as hhhp ~ hhph, 
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and phhh ~ hphh. The term phph is equivalent to three other terms: hphp, 
phhp and hpph. 

There are therefore only the following 9 possibilities, which are (excluding 
the factor ^^i^ which is the same for all of them) : 

(1) pppp =^ {i\ fej^p bs,p \k)Sj^i (2.31) 

(2) ppph ^ 2(-l)"^(^| blp b,,, \k){j\bs,,\l) 

(3) hppp ^ - 2(-l)"'=(^| b^p is,, \k){j\bl,\l) 

(4) pphh {i\ blp b^p \k){j\ b^,h hh 10 

(5) phph - 4(i| b^,p \k){j\ bl f^ bs,h 10 

(6) hhpp {i\ b^,p bs,p \k){j\ bl,^ b\,h 10 

{7) phhh 2{-ir^{t\blp\k){j\bl,L,,,bsM\l) 

(8) hhhp ^ - 2i-ir^{t\bsjk){j\ bl^ 6^^, b,,n 10 

(9) hhhh ^ 5i,k{j\ blh %,h Kh kh 10- 

We again remark that the enumeration over a, 7, 5 is done, in practice, 
only for cases (2) — (8), since for the cases (1) and (9) we use the operator 
which was already calculated as the required sum (two such operators exist, 
for both particle- and hole-states). 

Finally we arrive at the superblock Hamiltonian matrix. Its size can vary 
since we take only super-states with an equal number of particles and holes, 
and the number of such combinations depends on the states kept in each 
block, which can be changed in every iteration. The upper limit, of course, 
is 4p^, where p is the block size of both P and H. Anyway, we diagonalize 
the Hamiltonian, i.e., we find its ground state ^0 for which H\^q) — -EqI^o), 
and we move on to the next step. 

Step 4: The Density Matrix and Hilbert-Space Truncation 

Prom the ground state |^o) we build two density matrices, one for the P- 
block and the other for the H- block. Since actually |^o) encodes the indices 
for both blocks, i.e., 

1*0) = (2.32) 
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we can trace out half of these indices in order to build each density matrix. 
For the P density matrix we define (omitting the subscript from \l'o) 



The next step is to diagonalize these matrices and get, for each, the 
eigenvalues with the highest eigenvalues as a truncation operator. If we 
denote p\ua) — i^aWa) (for each one of the blocks), then we choose the 
eigenvectors \ua) with the highest cUaS to form an 2p x p matrix O. We then 
use the O matrix in order to transform all the enlarged operators (from step 
2) to their original small size p, as it was at the beginning of the current 
iteration. 

Now, of course, we are able to add another couple of states (one P and 
one H) to the game, and thus we return to step 2. 

As steps (2) — (4) are repeated iteratively, the number of states involved 
in the physical game increases, and the states we are considering are further 
away from the Fermi energy. One may thus expect that as this distance 
increases, the influence of the additional levels will decrease. 

A few minor remarks: in the first few iterations, when the block size is 
still smaller than p, the truncation is not done, the operators continue to the 
next iteration as they are, and the block size grows. When we are not at 
half fiUing occupation, the number of P-states is different from that of the 
H-states. For example when Up > Uh, after we finish adding all the H-states 
we continue in the process but add in each iteration only the next P-state. 
In the iterations that follow, the H-blocks will remain the same. 




(2.33) 



while p^'*) is defined as 




(2.34) 



Chapter 3 

A level coupled to a ID 
interacting reservoir 

In this chapter, we explore the ground state properties of a one- dimensional 
system, consisting of a gate-controlled dot coupled to an interacting reservoir, 
using the numerical density- matrix renormalization-group method. First, we 
calculate a few physical observables of the system in the thermodynamic 
limit and zero temperature in the non-interacting case. We then move to 
investigate the Tomonaga-Luttinger liquid phase of the interacting lead. We 
concentrate on the influence of interactions in the lead, as well as dot-lead 
interactions, on the width of the dot filling as a function of the chemical 
potential, and on the position of the dot level. 

We also study other phases of the lead, i.e., the charge density wave and 
the ferromagnetic phases. With particle-hole symmetry, it is known that for 
different values of the interaction strength the lead is described by different 
phases (see section [1.2. 3p . We show that a semi-infinite charge density wave 
coupled to the dot undergoes a quantum phase transition when the dot's level 
crosses the wire's chemical potential. On the other hand, in the ferromagnetic 
phase there is a simple level crossing at the same point. 

3.1 Introduction 

The properties of one- dimensional (ID) interacting systems have attracted 
much interest going back half a century ^Ol Hi], [38]. Much recent effort 
has concentrated on understanding the conductivity and I-V characteristics 
of a Tomonaga-Luttinger liquid (TLL) coupled to an impurity [l3j. These 
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properties are probed experimentally by measurements of the temperature 
dependent conduction through ID systems [Mj, and tunneling spectroscopy 
into ID wires [65J. The measurements of Ref. [65] also indicate a localization 
transition in the wire for low densities which might be associated with a 
charge density wave (CDW). Signatures of a CDW have been also observed 
recently, in finite ID wires coupled to dots, when strong magnetic fields were 
applied [55] . 

A generic model for the situation of a ID system which is coupled to an 
impurity, is a quantum dot (QD) coupled to a lead. If such a dot is controlled 
by a gate one can change its orbital energy, and measure physical quantities 
which are either on the lead or on the dot. For example, the occupation of the 
dot level is experimentally accessible by the charging effect of the impurity 
on a quantum point contact (QPC) in its vicinity [67] . 

3.1.1 Level Broadening vs. Conductivity 

The difference between measuring the conductivity through the dot-lead sys- 
tem or the local density of states at the impurity [68], and probing the dot 
occupation using, e.g., a QPC, must be emphasized. Essentially, as was noted 
in chapter [H any impurity which is placed in a TLL will lead to an insulating 
behavior [43]. The resonance conductance through a QD coupled to a pair 
of TLL leads (see Fig. 13.1b ) was found to produce infinitely sharp Coulomb 
blockade peaks at zero temperature [B]. Thus, no level broadening of the 
dot states is exhibited in the measurement of the conduction through that 
dot. 

Nevertheless, this does not imply that coupling a dot to a TLL has no 
effect on the width of the dot filling as a function of the chemical potential. 
Consider for example the arrangement depicted in Fig. 13.1b . A dot is con- 
nected to a TLL lead, while its occupation is measured by a QPC. Thus, the 
TLL acts as a reservoir for the dot, while the QPC is used to probe the dot's 
level broadening. In such an arrangement, any additional broadening of the 
levels due to the coupling to the reservoir will be seen in the shape of the 
conductance through the QPC. 

The conductance through the QPC in the geometry described in Fig. 13. lb 
is directly proportional to the occupation of the dot's orbital due to the 
capacitive coupling between the charge of the dot and the QPC (it is assumed 
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that no tunneling occurs between the dot and the QPC) [67] . Thus, in 
principal, the dot occupation ridot may be read off the conductance through 
the QPC and the effect of coupling of the dot to the interacting reservoir can 
be measured. 

The difference between the two arrangements depicted in Figs. 13.1b and 
13.1b is that while the first case (Fig. 13.1b ) essentially probes the enhancement 
of the backscattering in the TLL in the vicinity of the Fermi energy, the sec- 
ond (Fig. 13.1b ) explores the broadening of the level due to coupling to states 
which may be far from the Fermi energy. Therefore, one might expect the 
broadening of the level measured in the second arrangement to approach the 
conventional Breit-Wigner form, although some signature of the interactions 
is anticipated. 




Figure 3.1: Different coupling schemes of a QD and leads, (a) A QD 
coupled to two TLL leads represented by the wide lines, (b) A QD coupled 
electrostatically to a QPC through which the conductance in measured, 
and to an interacting reservoir, (c) The diagrammatic representation of 
the RPA approximation of the self energy. The line corresponds to the 
lead Green function, the black dot to the hopping into the dot and the 
wiggly line to the interaction. 



3.1.2 Chapter's Outline 

The rest of this chapter is organized as follows. In the next section, 13.21 we 
discuss the Hamiltonian model, and briefly describe the numerical methods 
we usel]. We then dedicate section 13.31 to an analytical calculation in the 
non-interacting case, in which we exactly formulate the physical properties, 
which we later calculate numerically for the interacting system, by using the 



detailed description is given in chapter [21 
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Green functions techniques. Next we present the numerical results for in- 
teracting systems. In section 13.41 we show the dot population in the most 
general TLL phase. We continue, in section [231 with a version of the Hamil- 
tonian which conserves particle-hole symmetry, and compare the system's 
properties between the TLL and the CDW phases. Three physical quantities 
are obtained: the dot occupation, the total system's population and the free 
energy. Focusing on the latter, we compare, in section [3^ its dependence on 
the gate voltage to that of the ferromagnetic (FM) phase, and show that in 
the CDW there is a quantum phase transition (QPT) when the dot's energy 
crosses the Fermi energy. In the FM phase, on the other hand, a simple 
level crossing occurs at that point. We continue the chapter in section [3171 by 
exploring another consequence of the particle-hole symmetry of the Hamilto- 
nian, which results in a non-trivial integral form. We finish by a summation 
of the main findings, and point out a possibility for a future research. 

3.2 Model 
3.2.1 Hamiltonian 

The system we investigate in this chapter is composed of a QD which is 
coupled to a ID lead. The electrons moving in the system are spinless, the 
QD is restricted to have a single level, and there are interactions between 
nearest-neighbor (NN) electrons in the lead. The Hamiltonian describing 
such a system can be written as H = Hdot + -^dot-icad + -^load, where ^dot 
and Hiend are given, as was explained in chapter [H by 

-f^dot = toa^a, (3.1) 

L-l L-1 

Aiad = -i + H.c.) + CjCjcj+iCj+i, 

and ^dot-lead, which connects the QD and the lead through hopping and 
dot-lead interaction terms, is 

-^dot-lead = -V{d^ci + c\a) + Id\a^ac\ci. (3.2) 

We denote the dot's energy level by eo, V (t) is the dot-lead (lead) hop- 
ping matrix element, and I^i (/) is the dot-lead (lead- lead) NN interaction 



3.2. MODEL 



55 



strength. Following our convention for the operators notations, (d) is the 
creation (annihilation) operator of an electron in the dot, and Cj [cj) is the 
creation (annihilation) operator of an electron at site j in the lead. The lead 
hopping matrix element, t, is taken as 1, in order to set the energy scale. 

In order to obtain a particle-hole symmetric version of the Hamiltonian, 
one should consider a positive background in the two interaction terms. 
Throughout this chapter, for the particle-hole symmetric version we will sub- 
stitute /di = in the _ffdot-iead term, and we will replace -ff/g^d 

L-l L-l ^ ^ 

ASd = + ^-^o + - 2^i^U^^+^ - 2)- (3.3) 

i=i i=i 

In the following we use both forms of these Hamiltonians. When the 
electron-electron interaction is not considered, as in section 13. 3[ the two 
forms -^/g^d ^^'^ ^leld identical. For the interacting case, however, they 
are different. In section 13.41 we examine the TLL phase by considering the 
non-symmetric Hamiltonian Hi = H^^t + -^dot-iead + Aiad- explore the 
interaction effect on the dot population in two cases. The first case is when 
an electrostatic potential between the dot and the lead does not exist, i.e., 
for Jdi = 0. In the second case we assume Jdi = /. In section [331 we compare 

between different phases of the wire by taking the particle-hole symmetric 

^ ^ ^ ^ (2) 

version of the Hamiltonian H2 = i/dot + -f^dot-iead + -f^iead ^i^h Jdi = 0. 

3.2.2 Diagonalization Method 

For each Hamiltonian model we discuss, the grand canonical potential Cl = 
H — fiNe, where Ne is the particle-number operator and /i is the chemical 
potential, is diagonalized using a finite-size density-matrix renormalization- 
group (DMRG) calculation [IHl [69] for different values of V^, J and eo or /i, 
and with a lead of up to L = 500 sites. For the cases in which a particle-hole 
symmetry is required, fi is set to zero. 

As was noted in chapter [2l the total number of particles in the system 
is not fixed during the DMRG process, so that the results describe the ex- 
perimental situation of a finite section of a ID wire which is coupled to a 
dot and to an external electron reservoir. In particular, the total occupation 
obtained for such a system can be non-integral. 
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The following ground state properties are calculated as a function of eo 
(or fi): the system's grand potential Q, the total number of electrons N and 
the dot population ndot- The population of lead sites are also calculated in 
order to differentiate between local effects of the dot population and global 
phenomena in the lead. 



3.3 Non-Interacting Case 
3.3.1 Exact Calculation 

In order to estimate the influence of the interactions in the lead on the 
different properties of the system, we should start from considering the non- 
interacting case. The coupling of the dot state to the continuum (akin to 
the Fano-Anderson model) may be treated using standard Green function 
technique [70] which leads to: 



1 Q'S(e) 

ndot(/^,eo) = - / 7 «?vr \\2^{c^Y^f \\2 ^^^ (3.4) 

J-oo (e - eo - 5tS(e))^ + (bS(e))^ 

where S(e) is the self energy given by: 

efc are the eigenvalues of the lead, Vk is the coupling between the eigenstates 
in the lead and the state in the dot and 5^0. 

For the idealized case, the density of states in the lead is constant (i.e., 
€-k = k/Lu, where u is the (constant) local density of states, and L is the 
leads length). The coupling is = ^fafLV (a is the distance between NNs), 
and under these conditions 

resulting in 3S(e) = T{av\y\^ = r/2 and 3ftS(e) = 0. Thus, one obtains the 
Breit-Wigner formula: 



2 



ridotifJ', eo) = - TT^de. (3.7) 

'-00 (e 
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For the tight-binding model discussed efc = — 2tcos(A;a) and = ^/2a/L sm{ka)V, 
resulting in 



S(e) = {V/tYe/2 + tiVyt)^/l-ie/2ty. (3.^ 

We thus find 



TT 



1 Vy^-w 

ridot(/x,eo) = - / — — -ode. (3.9) 



--¥(l-6) + [(l-S)^-^o]- 



The change in the total density of states in the system, due to the presence 
of the dot, can be found in a similar method to be 

Az/(e) = ln(e - cq - S(e)). (3.10) 

Therefore, the change in the number of electrons in the entire system is 

AN{n,ea)= r Au{e)de= (3.11) 

J~2t 

11 _i /i - ep - g?S(/i) 
2^^ ^^'^ ^SG^I) ' 

and the change in the free energy of the system is 

An(/i, eo) = - /iAA^(/i) = (3.12) 
r (e - /i)Az/(e)de = - AN{e,eo)de. 

3.3.2 Numerical Results 

The non-interacting case, in which the equations given above are relevant, 
gives an important possibility to examine the results of the numerical DMRG 
method. In Fig. 13.21 the numerical results of the DMRG for the dot popula- 
tion are presented in symbols, and are compared to the calculated formula 
Eq. fl3.9p . A very good correspondenc^^ is obtained for all values of the 
dot-lead coupling V. 



good correspondence is obtained for the other two formulas as well, i.e., between the 
DMRG results for the total number of electrons and for the free energy, to the prediction 



of Eqs. (j3.1ip and (|3.12p . These results are shown in section 1331 
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Figure 3.2: The dot population as a function of the chemical potential 
for eo = —1.1 is shown for various values of V, and without interactions. 
The symbols are the DMRG results and the lines are the exact formula, 
Eq. 

3.4 The Dot Population in the TLL Phase 

In this section we discuss the results obtained for the Hamiltonian Hi, con- 
sidering a constant eo and enumerating over /i. Alternatively, we will consider 
a constant /x, and enumerate over eo- In both cases the potential H — fiNe 
is diagonalized using the DMRG method for a lead of a few hundreds of 
sites, and the curves of the dot population, whether ndot(Ai) or ndot(eo); are 
calculated for several coupling strengths, V, and interaction strengths, I. 

3.4.1 Influence of Interactions in the Lead 

First we discuss the case in which Jdi = 0. Typical results for ndot(eo) and 
n^otifJ') are shown in Figs. l3.3 land l3.4[ As can be seen in both cases, changing 
the interaction strength in the lead (J) almost does not affect the level posi- 
tion, which is centered near the constant /i in Fig. I3.3[ and the constant eo 
in Fig. 13.41 The level width, however, decreases with increasing interactions. 
Nevertheless, it can be clearly seen that the shape of the curve is similar for 
all values of interaction strength shown. 

In order to understand the role played by the interactions in the lead, 
when it is described by the TLL theory, one must remember that the dot 
occupation ndot is determined by contributions from all energies, and the 
region around the Fermi energy, which has the special TLL behavior, does 
not play a unique role. Therefore, one could expect a simple perturbation 
description of the interactions in the lead to suffice. Indeed, the effect of the 
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Figure 3.3: (a) The dot population as a function of the level energy 
for V = 0.2 (dashed lines) and 0.4 (solid lines), and a constant chemical 
potential ^ = — 1. The interaction strength / takes values between and 
3.5, in jumps of 0.5. (b) The curves for 1 = 0, 3.5 are redrawn (symbols) 
together with the fit to the Eq. (j3.9p (lines) with an effective coupling V^s. 



electron electron interactions in the RPA approximation on the self energy 
(see Fig. 13.1b ) may be written as [7lJ: 

E(e) = xS°(e) (3.13) 

where S''(e) is the non-interacting self energy and 

Here we assumed a constant local density of states u (for the tight binding 
lead u = (a7rt)~^ which ignores local density of states variations), thus one 
obtains 5>S(e) = (VV[t(l + I /nt)])^! - {e/2ty and 3fJS(e) = {Vy[t{l + 
I/7it)])e/2, corresponding to replacing in Eqs.f l3.8p and (13. 9p by an "ef- 
fective" coupling 

Kff = V^V(l + /M)- (3.15) 

Returning to the results obtained by the numerical DMRG calculations, 
the curves of Figs. 13.31 and 13.41 can be now fitted to Eq. (13. 9p with V^g as 
a fitting parameter. It is easy to see that the Breit-Wigner form fits quite 
well even in the presence of strong interactions, and with the level position 
eo remaining constant. The effect of interactions is indeed limited here to 
a decrease of V^g, i.e., a decrease of the level width F. The values of V^^ 
extracted from the fit are plotted in Fig. 13.51 and compared with the RPA 
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-1.15 -1.1 -1.05 -1 

Figure 3.4: (a) The dot population as a function of the chemical potential 
for y = 0.15 (solid lines) and 0.3 (dashed lines), and a constant level energy 
eo = —1.1. The interaction strength / takes values between and 2, in 
jumps of 0.5. b) The curves for / = 0, 2 are redrawn (symbols) together 
with the fit to the Eq. ()3.9p (lines) with an effective coupling T4fj. 
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I 

Figure 3.5: Vcs/V'^ as a function of / (symbols) as obtained by fitting 
the ndot(eo) curves of V = 0.2,0.3,0.4 and 0.5 to Eq. The line 

corresponds to the dependence according to RPA prediction Eq. ()3.15p . 
and is in good accordance with all curves for small /. 

predictions of V^g = V^/i^l + I/nt). A rather good correspondence, especially 
for small values of /, is observed for all values of V. 

Thus, although the TLL has a vanishing local density of states at the 
end of the lead in the vicinity of the Fermi energy, the dot level is broad- 
ened, since all the reservoir states take part in the broadening mechanism. 
Nevertheless, as we have seen, the interactions in the reservoir influence the 
width of the resonance. One might gain some insight from the following 
consideration: For the non-interacting case (for constant density of states 
in the reservoir) the width is equal to F = 2TTai>\V\'^, which may be rewrit- 
ten as F = 2iT{a/ L)\V\'^dN/ dfi. The thermodynamic inverse compressibility 
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dN/ dfi is affected by the interactions [72] . Lets consider the compressibihty 
dfi/dN. In the lowest order approximation [71] dn/dN = {Lu)^^ + e^/C, 
where C is the capacitance of the system. For NN interaction = al and 
as usual C ^ L. Therefore, dN/d^ = Lv/{1 + avi). Inserting this to the 
expression for F, we get a result similar to the RPA approximation results in 
Eq. fl3.14p . Although capacitance is proportional to the length L of the lead, 
so is the density of states in a ID system, and therefore it has an influence 
even for an infinite lead. 

3.4.2 Influence of Dot-Lead Interactions 

In order to consider interactions between an electron occupying the dot, and 
the electrons in the lead, we now turn on the dot-lead interaction term, and 
we choose to take I^i = /• The corresponding n(yu) results (Fig. 13.61) clearly 
show a change in the resonance width, but also a change in the level position, 
which was absent in the previous case. Nevertheless, these results can still 
be fitted to Eq. (13. 9p . using also eo,efr as an additional fitting parameter, in 
addition to V^g. As can be seen, Eq. (13. 9p describes this system quite well. 
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^ -1.1 -1 -0.9 -0.8 

Figure 3.6: (a) Population of the dot and of the first site of the lead as 
a function of the chemical potential, for V = 0.15. Dot-lead interaction 
was included taking Idi = I. The curves shown are for / between and 3, 
in jumps of 0.25 (full lines) for the dot population, and for I = (dashed 
line) and 7 = 2 (dotted line) for the lead population, (b) The plots for 
I = 0, 1.5 and 3 (symbols) together with the best fit to Eq. (|3.9p . 

The movement of the resonance center, eo,eff, as well as the width, V^g, 
that were obtained from the fit are shown in Fig. 13.71 For small values 
of interaction both grow linearly. First lets try to explain the shift in the 
resonance center. As noted, almost no shift was seen for I^i = 0. In the 




62CHAPTER 3. A QD COUPLED TO A ID INTERACTING RESERVOIR 



presence of weak dot-lead interactions, one may approximate the dot-lead 
interaction operator by H^i ~ If^iriio^d, where rii represents the average 
occupation of the first site in the lead. As can be seen in Fig. 13.6b . ni is 
not very sensitive to the occupation of the dot and may be replaced by its 
typical value. Thus, the energy of the orbital in Eq. ( 13. ip may be rewritten 
as eo,efr = Co + ^i-^di- Indeed, this formula fits well the numerical results for 
small values of Idi (as can be seen in Fig. 13.7( a)). for rii = 0.14. This result 
agrees well with the value rii ^ 0.15 in the region of the resonance, taken 
from the data of Fig. 13.6( a). 



Figure 3.7: (a) eo,eff and (b) Vjj as functions of Idi = I (symbols) and 
linear fits for the region / < 0.5 (lines). 

A more striking feature is the behavior of V^g, i.e, the width of the reso- 
nance. There is a distinct qualitative change in the width behavior, compared 
to the case without dot-lead interactions. As opposed to the monotonic de- 
crease of Vgg, which was demonstrated in Fig. 13. 5[ Fig. 13.7b (symbols) shows 
that V;| increases with J, until a maximal value is achieved around 1 = 2. 
For larger values of interaction a decrease in the width is observed. 

This enhancement of V^g is associated to the interplay between the popu- 
lation of the dot level to the depopulation of the first site in the lead, ignored 
in our treatment of eo,efr- This leads to a reduction in the effect of the dot-lead 
interaction which results in an increase in the width as depicted in Fig. 13.7b . 
For weak interactions the enhancement of V3f is linear. 
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3.5 A Comparison Between the TLL and the 
CDW Phases 

The ID wire can be, as noted in chapter [H in other phases than the TLL, 
e.g. the CDW phase. The cleanest appearance of the CDW phase occurs 
at /i = 0, when the system is half-filled and particle- hole symmetry exists. 
Therefore, in order to compare the results for these two phases, we use, here, 
the symmetric version of the Hamiltonian, H2. 

In order not to increase the number of free parameters in this section, the 
model discussed here does not contain a dot-lead interaction term. Neverthe- 
less, as the results of the previous section have pointed out, and according to 
the results we are going to present, we note that such a term is not expected 
to influence qualitatively our results. 

3.5.1 Level Occupation 

We begin by comparing the behavior of ridot in the two phases of the inter- 
acting lead. As pointed out in section 11.21 the system experiences a phase 
transition (of a Kosterlitz-Thouless type [32]), between the TLL and the 
CDW phases, at / = 2t. Indeed, the results presented in Fig. 13.81 show a 
qualitative difference between ndot(eo) in these two regimes. We have shown 
in the previous section that in the TLL phase the curves fit quite well the 
Breit-Wigner formula with an effective coupling Veg. In the CDW phase, 
however, the width of the level becomes zero, and a jump in n^ot occurs at 
69 = 0. This jump is associated with the degeneracy of the CDW ground 
state. 

In order to support this argument, we first calculate the properties of the 
ground state for a lead not coupled to a dot, by taking V = t and the limit 
Co ^ 0. The CDW order parameter can be defined as [73] 

P{i) = i(_i)^[2n(i) - n{i - 1) - n{i + 1)], (3.16) 

where n{i) = {c\ci) is the occupation of the i'th lead site in the ground 
state. The value of P{i) does not change much as a function of the spatial 
coordinate, except in the vicinity of the lead's edges. We thus define P = 
P{L/2). 
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Figure 3.8: The dot population ndot as a function of the level energy 
eo- For I < 2t (TLL phase) the curves fit the non- interacting formula with 
an effective coupling constant Vcs (lines - DMRG results, symbols - fit to 
Eq. (|3.9p ). while for I > 2t (CDW phase) the width is zero. Exactly at 
the transition point (/ = 2t) n^ot is continuous but doesn't fit the non- 
interacting formula. 



In Fig. I3.9l the dependence of P on the interaction strength J, as obtained 
by the DMRG method, is shown, and compared to the exact results [35] . For 
/ < 2t, the system is in the TLL phase and indeed P = ± 10~^. In this 
case the population of each lead site is 1/2 and the lead is half filled; i.e., 
N = L/2. 

As expected, the metal-insulator transition occurs at / = 2t. For / > 2t, 
in which the system is in the CDW phase, the value of P is finite. For values 
of / which are far from the transition point (/ > 2.5t), we get a very good 
fit to the theory. Near the transition point the numerics tend to emphasize 
the charge oscillations, resulting in too large values of \P\. This tendency, 
however, does not affect the following qualitative conclusions. 

As can be seen, the CDW ground state is different for the cases eg 
and eo — > 0~, resulting in two values of P (i.e., ±|-P|) for any value of / in 
the CDW regime. For eo = the ground state is two-fold degenerate in the 
thermodynamic limit. 

In the CDW case, special care should be devoted to the number of elec- 
trons in the system, and to the difference between even and odd lead length 
L. For eQ ^ and V = t, the lead length is effectively L + 1. Denoting the 
ground state for eo 0^ by and its population in each site as n^j^\ 

where j = for the dot and I < j < L foT the lead sites, the particle-hole 
symmetry implies that n^j ■* = 1 — rij^^ everywhere. For odd L (even L + 1) 
one has the requirement -* = n^i^j, and the number of electrons in both 
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Figure 3.9: The CDW order parameter P as a function of the interaction 
I for a 300-sites lead (DMRG - symbols, theory - lines). For the CDW 
phase P is inverted between the cases eo — >■ 0"*" and €q ^ 0~, while in the 
TLL phase both cases result in P = 0. Inset: P as a function of eo for 
/ = 1.5 (triangles), 2.5 (squares) and 3.5 (circles). The TLL case results 
in a constant P = 0, while an inversion of P occurs at eo = for the two 
CDW cases. 



and is = (L + l)/2. When L is even (L + 1 is odd), however, 
one has an additional symmetry requirement, nj = nL_j, for both and 
'ipQ~\ and the states differ in their electrons number: i/jq'^^ contains N = L/2 
electrons, for any / > 2, while ipQ ^ has iV = L/2 + l. It is worth noting that 
nS+^ > 1/2, while n^"^ < 1/2. 

When eo ^ 0, i.e., a dot with a finite on-site energy is connected to the 
lead, it has some local infiuencf f] on the ground state in its vicinity. For a 
lead in the CDW phase, the two states described above are slightly modified, 
and can be denoted by tjj^^^ and Nevertheless, the main influence of 

the dot is to lift the degeneracy between those states. If eo < 0, the dot state 
population is high and is energetically preferable due to the dot-lead 
hopping term. For eo > the opposite happens, resulting in a preference 
of ip^~^\ The occupation of the dot-lead system in the CDW phase is shown 
in Fig. 13.101 When the dot's orbital energy changes from a negative to 
a positive value, the system switches from ip^~'^ to ip^~^\ In the following 
section we show, by a calculation of the free energy, that this is indeed a first 
order QPT. The resulting phase diagram is drawn in the inset of Fig. 13. Ill 

3.5.2 Total Number of Electrons 

In Fig. 13.111 the total number of electrons in the system is presented. For 
I = 0, N fits the predicted formula, Eq. (13.111) . well. As noted above. 



^The investigation of this influence takes a major part of chapter |4l 
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Figure 3.10: The occupation of the dot and the first lead sites of the 
CDW states tp^^^ (circles) and tpQ ^ (squares), of a lead with 300 sites and 
7 = 4 which is coupled to the dot (filled symbols) with V = 0.8. For 
eo = these two states are degenerate in the thermodynamic limit. Lines 
are guide to the eye. 



the agreement between the numerical results for a finite lead and the exact 
result for a semi-infinite lead, is obtained due to the fact that the number of 
particles can vary during the DMRG process. For the TLL phase (J < 2t), 
A^(eo) looks quite similar, varying between L/2 + 1 (at eo — oo) and L/2 
(at eo oo), taking the average value (L + l)/2 at eo = 0. As in the case of 
ndot, the results fit Eq. f l3.1ip with a renormalized dot-lead coupling for small 
values of I. Increasing / towards the transition point (/ > 1.5t), results in a 
less accurate fit. 




Figure 3.11: The total number of electrons in the system, as a function 
of eo, for L = 299, n = 0, V = 0.3, and different values of /. Symbols 
- DMRG results. Lines - fit to Eq. ()3.1ip . Inset: the phase diagram: a 
Kosterlitz-Thouless transition (solid line) occurs at / = 2, and a first order 
phase transition (dashed line) at eo = for / > 2. 



The values of Vcs{I) obtained by the fit of the curves to Eq. (13. lip , are 
in good agreement with the values obtained by the fit of ridot to Eq. (13.91) . 
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We find that Ves decreases monotonically with increasing J, exhibiting the 
RPA hke behavior described in section I3.4[ 

The CDW phase (/ > 2t), however, is quahtatively different: for an odd 
lead length A^(eo) remains a constant integer (A^ = (L + l)/2) which does 
not depend on eo at all. For even L, N = L/2 + 1 for eg < 0, and = L/2 
for Co > 0. This is a direct result from the switch of the system ground state 
from -ip^'^ to 'ip^^^ at eo = 0. For an odd lead length, the total number of sites 
in the system (L + 1) is even, so that is equal for both states. For an even 
lead length, L + 1 is odd, and the total number of electrons is changed by 
one when eo passes 0. Except for the decrease of one electron at eo = for an 
even lead, A^ remains constant in the CDW phase, independent of eo. Thus 
even with the continuous change in the population of the dot as a function 
of the dot's level for eo 7^ 0, the number of electrons in the entire system 
remains constant. The change in the occupation of the dot as a function of 
eo is compensated by the lead. 

This difference in the behavior of A^ as a function of eo between the two 
phases (i.e., constant for the CDW phase, compared to a continuous decrease 
for the TLL phase) is a direct manifestation of their transport properties. 
The TLL phase is metallic, and therefore compressible. Hence, infinitesimal 
changes of the electrons number are possible. On the other hand, the CDW 
phase is insulating and thus incompressible, which results in a constant A^. 

3.5.3 Free Energy 

In Fig. 13. 12^ typical numerical results for the free energy ^(eo) — ^^(0) as a 
function of eo are shown. A perfect fit between the DMRG results and the 
exact formula Eq. (13.121) for / = is obtained. In the TLL phase (/ < 2t) our 
DMRG calculations show that the effect of interactions on Q can be fitted 
by replacing V in Eq. (I3.12p by the same effective coupling Ves obtained for 
the behavior of Udot and of A^ which were discussed above. For the CDW 
phase (/ > 2t), however, there is obviously a qualitative change in the energy 
curve: the dependence of Q on eo is linear both below and above eo = 0, with 
an abrupt change of ^ at eo = 0. 

These results point out that the single impurity, connected at one end of 
a long interacting lead which is in a CDW state, has a well defined influence 
on the ground state of the entire coupled system. As discussed above, when 
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Figure 3.12: The ground state free energy as a function of eo for different 
interaction strengths. The lead-dot couphng was taken asV = 0.3, and the 
interaction strength / takes values between and 3.5, in jumps of 0.5. The 
lines represent the DMRG results (dashed line -1 = 2), while the symbols 
show the fit to Eq. (j3.12p for values of / lower than 2. Inset: ri(eo) — cq) 
as a function of eo for / = 0,1,2,3 (symbols) and f2(eo) — il(— eo) = eo 
(line). For all values of interaction we get |l^(eo) — i7(— eo) — eo| < 10~^. 
This result is discussed in section 13.71 

the dot level passes through eo = 0, the lead's population is inverted at 
every site, leading to an inversion of the CDW order parameter, as presented 
in Fig. I3.9( inset). The dot population is inverted as well. As a result, a 
dramatic change in the dependence of the free energy of the system on eo 
occurs. 

3.6 Quantum Phase Transition 

3.6.1 QPT or a Simple Level Crossing? 

Since Q is the free energy of the system, the jump in its first derivative might 
be a sign of a first order QPT. In order to see whether this non-analyticity 
of the free energy is not just a trivial level crossing (LC) of two levels in 
the system, the dependence of the transition shape on the system size L is 
explored. For the case that the sharp transition in the energy results from the 
fact that the external field (in our case the gate voltage) commutes with the 
Hamiltonian, and thus a LC is possible, no size dependence of the sharpness 
of the transition is expected. On the other hand, a real QPT will become 
sharp only in the thermodynamic limit (i.e., semi-infinite lead). 

In order to compare the two scenarios (LC vs. QPT) we solve the same 
Hamiltonian with strong attractive interactions, i.e., / < —2t. As discussed 
in chapter [H such strong attractive interactions regime corresponds to the 
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FM phase, while the CDW is equivalent to the Neel phase. In the fermionic 
case the FM phase yields a two-fold degenerate ground state, composed of 
states in which the sites are either entirely occupied or emptjQ. Pinning by 
the dot lifts the degeneracy, because an occupied dot causes a preference of an 
empty lead, and vice versa. Thus the influence of the dot on the coupled wire 
is superficially similar to the CDW case in so far as in both cases the ground 
state of the entire system is determined by the dot orbital. Nevertheless, the 
FM system is clearly a LC system, and no dependence on the system size is 
expected. 



a 
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Figure 3.13: The ground state free energy as a function of eo for different 
CDW system sizes, with V = 0.3 and I = 2.5t (dotted, dashed-dotted and 
dashed hnes). Inset: the results for L = 29 are compared to the case 
of / = —2.5t, which is inside the FM phase, and to the CDW case with 
L = 30. 



The dependence of f2 on eo was calculated for L = 29, 99, 299 in the CDW 
case, and for L = 29 in the FM case. From the results shown in Fig. 13.131 
(and its inset) it is clear that these two cases are different. The FM system 
is indeed a trivial LC system, with two competing states whose energies 
cross each other for eo = 0. Since these two states are eigenstates of the 
Hamiltonian even for a finite L, the size does not play a role, so that even 
for L = 29 one can see a sharp transition between the two ground states. 

On the other hand, for a small (but with an odd L) CDW system Q shows 
a smooth dependence on eg, without any non-analyticity. As a matter of fact, 
for any finite system the electron levels are expected to be mixed, resulting in 

"^No explicit restriction to a certain filling factor is assumed. 
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avoided crossing of the two lowest many-body levels [71] . In other words, the 
CDW states are not true eigenstates of the Hamiltonian for a finite system. 
Indeed, indications for the sudden jump of ^ are seen only for larger system 
sizes (i.e., for L > 200). 

The dot-dashed line in Fig. I3.13( inset) shows that for a short lead of an 
even size {L = 30 in that case) there is, however, a non-analyticity in f2 as a 
function of eo. The comparison between even and odd lead sizes emphasizes 
our conclusion stated above. While for an odd lead size ip^~'' and ■ip^'^^ have 
the same number of electrons, for an even size of the lead, the transition 
from to tjj^'^^ at eo = involves a decrease of the electrons number by 
one. The Hamiltonian H conserves the number of particles, so that in the 
case of L even, these two states are not coupled, and the transition between 
them is a simple LC, not showing a size dependence. For an odd L, however, 
and ip^^"^ are coupled by H, so that for a finite L they are actually 
mixtures of the CDW states, thus presenting sharper fi(eo) dependence for 
larger systems. We thus conclude, based on these two comparisons, that in 
the case of a CDW with an odd L this transition is a QPT, which happens 
for L oo, i.e., for a semi-infinite lead. 

3.6.2 Scaling Results of Small Systems 

As pointed out above, in the numerical results for large lead sizes, the non- 
analyticity of the free energy is seen clearly, so that it is obviously a first order 
transition. However, a similar conclusion can be drawn from the results of 
short leads, by scaling the results with the lead size L. Although for finite size 
systems the order parameter changes continuously, the slope of this change 
grows with the system size L. For a first order transition in d dimensions 
one expects a power law dependence of the slope as L'^, while for a second 
order transition the power law should be fixed by some universal exponents 

For the CDW model, the population of the dot plays the role of an order 
parameter, and one can check its behavior near eo = for different lead sizes. 
In Fig. 13.141 the occupation as a function of eo is shown for some short lead 
sizes {L = 29,49,69) and for a large one {L = 299). It is indeed seen that 
while for a long enough lead there is a jump in the dot population, for short 
leads there is a continuous change in the occupation, with a larger slope for 
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larger lengths. Zooming into the regime of eo = we find that although 
L > 200 is required in order to see a clear discontinuity in Q, L > 100 is long 
enough for showing a jump in ridot- The dependence of the slope on L, for 
L < 100, is shown in the inset of Fig. I3.14[ As one expects from finite size 
scaling predictions for first order phase transition, there is a very good linear 
fit of the the order parameter on the lead size. 




'0 



Figure 3.14: The dependence of the dot occupation on eo for different 
CDW system sizes, with V = 0.3 and / = 2.5. Inset: the absolute value of 
the slope in the limit eo — > as a function of L for short leads, together 
with a linear fit. 



3.7 Another Implication of Particle-Hole Sym- 
metry 

We now turn to another consequence of the particle-hole symmetry. For 
Co 7^ the symmetry results in the fact that the ground states for +eo 
and for — eo and ip^~^ respectively) have an inverted population ev- 

erywhere. It is thus clear that {ip^'^^Hiea.dli'^^^) = {ip^~^\Hies,d\'4'^~^) , 
similarly (?/'*^+) |-f^dot-iead|V'*-"^^) = (V^^"^|-ffdot-iead|^^~'')- This implies 

fi(eo) - ^^(-eo) = (3.17) 

eo<ot - (-eo)%ot = eo, 

where n'^^^^ ('^dot) represents the dot populations of ip^~^^ The last 

equality results from the symmetric population n^^^ = 1 — n^^^. 
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The relation Eq. f l3.17p does not depend on the interaction, and should 
exist for both the TLL and the CDW phases. Obviously it should be obeyed 
for the non- interacting case, where the energy is given by Eq. (13.121) . thus 
leading to the following nontrivial relation 



in which the result of the integral on the LHS does not depend on the pa- 
rameters V and t. 

A physical insight into Eq. (13.181) can be gained by taking the derivative 
of both sides with respect to eo, and using the definition of the self energy 
given above in Eq. (13.81) . This yields 



which is evident since the LHS is the occupation of the dot when jj > 2t. 

In Fig. I3.12( inset) a plot of fi(eo) — f2(— cq) as a function of eo is shown, 
for values of / between and 3. As can be seen, Eq. (13.171) is valid for all 
values of /. 

3.8 Conclusions and Future Prospects 

In conclusion, we have shown in this chapter that the occupation of the dot 
level can be used to identify the different phases in the wire. When the 
reservoir is in the TLL phase, the occupation of the QD does not show a 
jump when the dot level crosses the Fermi energy. The interactions in the 
reservoir, in that case, leave clear fingerprints on the width of the resonance. 
The main influence is a decrease in the resonance width due to a change in 
the inverse compressibility of the reservoir. On the other hand, a dot-lead 
interaction shifts the resonance position and may also enhance the width. 

In the two other phases, the CDW and the FM phases, the dot level splits 
the ground state degeneracy by favoring one of the two states depending on 
whether the QD level is empty or filled. Nevertheless, we have shown that 
the physics in these cases is different. While for the FM phase a simple LC 
occurs, with a sharp jump in the occupation of the level for any length of the 
wire, in the CDW phase the position of the dot level drives a first order QPT 




(3.18) 




(3.19) 
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in the thermodynamical limit between the two CDW states, while for a finite 
wire the jump is smeared. This phase transition shows all the hallmarks of 
a first order QPT, such as a size dependence, a jump in the order parameter 
and a discontinuity of the derivative of the grand canonical potential. 

* * * 

One of the results discussed in this chapter is the abrupt jump of Udot 
in the FM and the CDW phases, which results from the degeneracy of the 
ground state at eo = 0. In the TLL phase, as noted, no such a jump was 
demonstrated, and the dot occupation has a smooth analytic form. Nev- 
ertheless, Furusaki and Matveev have predicted [76] that a jump in the dot 
population should appear in the TLL phase as well, if the parameter k, which 
describes the interaction in the TLL (see section [L3l) . is smaller than |. 

In the model for a ID lead which we treat in this chapter, that of spinless 
fermions with a NNs interaction term, the TLL phase is restricted to the 
regime — 2t < I < 2t, and the relation between the interaction strength / 
and the TLL interaction parameter k, is given, for a half filled lead, by [T^ 
1^ = n — -i7 r/^o^M ■ Therefore, the values of k are limited to the regime k > h 

2 cos '-[—I/(2t)\ ' <=' 2 

At the point in which k approaches | the lead forms a CDW, and is no longer 
described by the TLL theory. Therefore, by using the NN interactions model 
of spinless fermions, with a half filled lead, one cannot get the Furusaki- 
Matveev jump. 

When the lead is not half filled, however, a CDW phase is usually not 
formed, even for strong interactions, since such a transition is based on umk- 
lapp processes which occur only for commensurate fillings. Whereas there is 
no analytic relation between k and / for systems which are not half filled, 
yet numerical studies of Haldane [7S] point out that a regime in which n < ^ 
exists, when the interactions are strong enough, and the lead occupation is 
not too far from half filling. 

Using a numerical method, however, it is a bit difficult to obtain this 
regime. The ground state of the TLL phase, with strong interactions (/ > 
2t) and near half filling, is basically built upon an half filled CDW, with 
some extra electrons which convert its behavior to a TLL. The number of 
these extra electrons is very important, since the TLL description assumes 
a linearization of the energy curve near the Fermi points (see chapter [1]), 
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and for a small number of particles it is not a good enough approximation. 
Therefore, unless one uses a very long lead (which might be computationally 
difficult), the strong interactions can cause the CDW ordering to be more 
favorable than a non-accurate description of the lead as a finite TLL. As a 
result, the TLL description, which is true for an infinite lead, is replaced by 
a finite-size CDW. 

There are also some modifications of the ID lead model which may be used 
to explore the n < ^ regime. For example, one can enlarge the interaction 
range, e.g. by considering also next-nearest-neighbor (NNN) interactions. 
Such a model was used in some previous studies, especially at half filling 
[79| [80l [81] , in which it results in three different phases as functions of Vi 
and V2 (denoting the NN and the NNN interaction strengths, accordingly). 
For Vi ^ V2 the ground state is a two-fold degenerate 2kF CDW, i.e., the 
occupation seems like ■ ■ ■ • o • o ■ ■ ■ . The other limit, in which V2 ^ Vi, 
results in a four- fold degenerate kp CDW creation, with sites occupancy such 
as ■■■••oo«»oo---. The intermediate regime, including the frustration 
line Vi = 2V2, is described by a TLL, at least for not very strong interactions 

m- 

As in the NN model, the half filled case in the NNN model might also 
be limited to k > |, as the numerical studies suggest. Nevertheless, when 
the system is not half filled, the case of k < | can be reached by this model 
with much smaller values of interaction [82j, in which the CDW excitations 
are weaker, thus shorter leads might suffice. 

Another possibility, which seems the most promising, is to consider spin 
1/2 electrons. In this case the interactions parameter has the analytic form 
^spmi/2 _ ^spmiess^2^ resulting, for a half filled lead, in the regime k > j. 
With this model any non-vanishing value of repulsive interactions result in 
K < |, and can demonstrate the Furusaki-Matveev jump. 

The required modifications of the DMRG process for NNN model or for 
spin 1/2 electrons are discussed in section 12.4.31 Here we just note that 
whereas in the spinless NN model there is an additional one state in each 
iteration, these modifications demand an iterative addition of a couple of 
states, whether they are two sites with spinless electrons, or a single site 
which might be doubly occupied. Therefore, in order to obtain the same 
accuracy, the DMRG block sizes, and thus the computational resources which 
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are required, are larger. The search for the Furusaki-Matveev jump is thus 
left for a future study. 



Chapter 4 

A level coupled to a disordered 
ID interacting reservoir 

In this chapter wc calculate numerically the Friedel oscillations caused due 
to an impurity located at one edge of a disordered interacting quantum wire. 
The electron density in the system's ground state is determined using the 
density-matrix renormalization-group method, and the Friedel oscillations 
data is extracted using the density difference between the case in which the 
wire is coupled to an impurity and the case where the impurity is uncoupled. 

For the non-interacting case, we develop an exact formula of the Friedel 
oscillations in the one-dimensional tight-binding model. The excellent fit to 
the numerical results serves as a proof for the accuracy of the method. 

For a one-dimensional wire which is described by Tomonaga-Luttinger 
liquid theory the oscillations of a clean interacting sample decay as a power 
law. We show that once the wire moves into an Anderson insulator phase, 
which happens due to the introduction of a disorder, the power law decay is 
multiplied by an exponential decay term due to the disorder. Scaling of the 
average Friedel oscillations by this exponential term collapses the disordered 
samples data on the clean results. The decay length is shown to decrease as 
a function of the interaction strength. 

However, when a short enough mesoscopic wire is in a charge density wave 
phase, the presence of weak disorder may not destroy the long range order, so 
that the wire will be described by a disordered Mott insulator phase. In this 
case we find that the effect of interactions is the opposite, and the disorder 
significance decreases as the interactions strength increases. 
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We prove that the length scale governing the exponential decay, in the 
Anderson insulator phase, may be associated with the Anderson localization 
length and thus be used as a convenient way to determine the dependence of 
the localization length on disorder and interactions. Our results, which show 
a decrease of the localization length as a function of the interaction strength, 
are in accordance with previous predictions. 

4.1 Introduction 

The interplay between repulsive interactions and disorder in low- dimensional 
systems, and their influence on the conductivity, were the subjects of many 
studies in recent years. Some of this interest was motivated by the exper- 
imental observations of a crossover in the temperature dependence of the 
conductance of low density two-dimensional electrons from an insulating- 
like dependence at low densities to a metallic one at higher densities |83j . 
Nowadays it is generally accepted that even if such a two-dimensional metal- 
insulator transition exists, it must be related to the spin degree of freedom 
jM] and therefore absent for spinless electrons. 

It seems therefore clear that for spinless one-dimensional (ID) systems no 
metal-insulator transition is expected for repulsive interactions, although for 
a certain range of attractive interactions a delocalized regime was found in 
several studies |85] . Nevertheless, it was shown that there might be a certain 
strong disorder and interaction regime, in which there is an increase of the 
localization length ^, defined through the dependence of the zero tempera- 
ture conductance on the system size: g{L) ~ exp(— L/^). Similar behavior 
was demonstrated also in properties which are usually related to ^ such as 
the persistent current [86]. A sample dependent increase in the localization 
length was also reported for weaker values of disorder and repulsive interac- 
tions for longer (of order of 100 sites) wires [S7j . 

On the other hand, several analytic studies which were performed in the 
Tomonaga-Luttinger liquid (TLL) framework [31] have concluded that the 
localization length of a ID wire decreases monotonically with increasing re- 
pulsive interaction. Using either renormalization group [88] or self consistent 
Hartree-Fock [HS] methods it was shown that the localization length, renor- 
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malized by the interaction, scales as 

^k) ~ {^oY^^'-'-\ (4.1) 

where is the locahzation length of the free electron system, and k is the 
TLL interactions parameter (see chapter [1]) with k = 1 for non-interacting 
electrons. Since for repulsive interactions k decreases as a function of the 
interaction strength, one finds that the localization length always decreases 
as a function of the interaction strength. 

One must be careful though to differentiate between weak and strong 
interaction strength. A clean ID system of spinless fermions on a lattice 
undergoes a metal-insulator phase transition between a TLL and a Mott 
Insulator (MI) as a function of the interaction strength. The MI phase, 
for strong interactions, appears for spinless ID electrons as a 2kF charge 
density wave (CDW). This phase transition, caused by umklapp processes, 
is exhibited for commensurate fillings. Once disorder is turned on, the TLL 
transport properties change drastically. For more than a decade it is well 
known [13] that the conductivity of a TLL wire vanishes in the presence 
of impurities, thus a metal-insulator transition as a function of interaction 
strength no longer exists. 

Nevertheless, a difference between the TLL and the CDW phases may 
still exist, even in the presence of disorder. When the interactions are weak, 
so that the wire is in a TLL phase, the addition of disorder turns the metallic 
system into an Anderson insulator (AI). However, for strong interactions (i.e., 
when the clean system is a CDW, which is a MI) the exact effect of disorder 
depends on its strength, and in general is not completely understood. While 
for clean systems the MI phase is a well studied problem [90] , the addition of 
disorder opens a few questions, which have attracted several studies in the 
last decade [HB ESI ESI El ES] . 

When the disorder is strong, i.e., when the random potential felt by the 
electron is much larger than any other energy scale in the problem, the MI 
state is destroyed, and an AI phase emerges. For weak disorder, however, it 
was shown in several studies that the Mott energy gap vanishes only when a 
finite disorder is introduced, so that below this critical disorder the MI phase 
is stable [^Sl EZ] . Usually this is not the case for a MI consisting of spinless 
particles, such as the CDW we study, since an Imry-Ma type of argument 
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[98] shows that the long range order is destroyed even for an infinitesimal 
disorder [9^. Yet, for a finite sized mesoscopic sample, the Imry-Ma length 
scale might be a few orders of magnitude larger than the sample's size, so 
that the effective ground state for a weak enough disorder remains a MI one. 
Increasing the disorder above a critical strength changes the MI state either 
to a Mott glasfl or to an AI [MI95] . 

Therefore, a finite size CDW state is expected to remain stable against 
the application of a weak enough disorder, i.e., to remain a MI state. For 
example, previous numerical simulations have presented the long range order 
of such a weakly disordered CDW [91j. Furthermore, such finite ID wires 
coupled to dots have been recently manufactured, and signatures of a CDW in 
strong magnetic fields have been observed f66] . In order to verify numerically 
the existence of a CDW order in the presence of disorder for the length scales 
considered, one can check the electron density of the entire system. 

Most studies on disordered ID wires concentrate on either the AI or the 
MI phases, thus a full comparison between the two regimes is still lacking. 
Nevertheless, a qualitatively different behavior between these two regimes 
was demonstrated in a few cases. For example, the effect of interactions 
on the persistent current in ID disordered rings was calculated in previous 
works [HU llOOj . and an important difference between the AI and MI phases 
was found. While for strong interactions and weak disorder (MI phase) the 
persistent current was reduced, for strong disorder (AI phase) an increase of 
the current was found. However, the exact diagonalization techniques which 
were used in these studies, are applicable only for very small system sizes. 

In this chapter we investigate the influence of interactions on the Friedel 
oscillations (FO) in a disordered ID wire, and compare this behavior between 
the AI and the MI regimes. We study interacting spinless electrons conflned 
to a ID wire which can be in either its AI or MI phases. The entire parameter 
range is very relevant to the explanation of recent transport measurements 
in various ID or quasi ID systems, such as single-wall carbon nanotubes 
[23l|2H[25], multi-wall carbon nanotubes [27j, polymer nanoflbers [28] and 
MoSe nanowires [30], which have led to many theoretical works [lOli 11021 
I103[ 11041 1105] . With this in mind, we study the effect of the interplay of 

^This issue will be addressed at the last section of this chapter. 
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interactions and disorder on the behavior of the FO in a wire due to its 
couphng to an impurity at its edge. 

4.1.1 Friedel Oscillations 

Once a single-level impurity (dot) is coupled to a clean metallic system the 
density of electrons in its vicinity oscillates with a 2kF period, and the en- 
velope of the oscillations decays as a power law of r, the distance from the 
impurity |106j . For non- interacting systems the perturbation of the density in 
the vicinity of the impurity depends on the dimensionality, d, of the system, 
and can be expressed as 

Hr) = aS2!P^£^, (4.2) 

where the coefficient A and the phase shift t] do not depend on r. These 
oscillations are the famous Friedel oscillations, which have been observed ex- 
perimentally during the last decade using various techniques, such as scan- 
ning tunneling microscopy in low temperatures [107] and X-Ray diffraction 

Whereas for higher dimensions {d > 2) Eq. (14 .2^ is in general true even 
in the presence of interactions, this is not the case for ID systems. For the 
TLL phase, using field theoretical approaches, it was shown [109j that the 
dependence is replaced by a different power law, x~'^. For the non-interacting 
case K = 1, it leads to the expected decay, while for repulsive interactions 
K < 1 and thus a slower decay of the FO envelope is expected. 

Therefore, the density change in the TLL phase, which is metallic, shows 
FO with a 2kF wave vector and a power law decay. In the CDW phase the 
picture is different, since the CDW phase is insulating, and thus the power 
law is replaced by an exponential decay. The length scale of this exponential 
decay is related to the CDW correlation length [32], C- 

In the presence of a weak disorder, one can expect to find a similar decay 
to that of the clean case, with an additional exponential decay due to the 
disorder. By calculating this exponential decay we are able to present a clear 
picture of the dependence of the decay length due to disorder on interac- 
tions, in both the AI and MI regimes. In the following sections we show a 
different behavior of this decay length in the two phases, and the origin of 
this difference is explained. 
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From Eq. (14.21) it is clear that the observation of the density fluctuations, 
either experimentally or numerically, is easier in the vicinity of the impurity. 
When disorder is also introduced, this distance becomes even shorter since 
there are also density fluctuations caused by the disorder. Yet, in common 
experimental ID situations disorder is usually present. Therefore, although 
the presence of disorder hampers observing the FO, it is beneficial to develop 
a method to tease the FO out of the density fluctuations of a disordered 
system. 

4.1.2 Localization Length in the TLL Phase 

In general, the study of a dot (or impurity) coupled to a ID lead, has been 
shown to shed light over the physics of the lead. Certain thermodynamic 
observables, such as the occupation of the impurity level (see for example 
chapter [3] and Ref. |76]) and the corresponding electron density changes in 
the lead 1111] , were recently used to analyze different wire properties, 
such as the strength and form of the interactions, and even the wire's phase 
(e.g., TLL vs. CDW). In a similar fashion, we show in this chapter how 
the electron density of a disordered wire in the TLL phase, coupled to an 
impurity level, can be used in order to extract its localization length. 

In the TLL case, for weak enough forward scattering we will show in the 
following that the decay length of the density- density correlations is equiv- 
alent to the localization length. It is important to note that the extraction 
of the localization length for interacting systems is plagued with difficul- 
ties. The straightforward method of measuring the decay length of the enve- 
lope of the single-electron state has no direct translation to a many-electron 
state. Nevertheless, one would prefer to stick to a ground-state property of 
the system, since the calculation of excited state dependent properties such 
as the conductance is computationally taxing. The sensitivity to bound- 
ary conditions (i.e., persistent current) which is the natural candidate for a 
ground-state property is problematic since it incorporates both interaction 
corrections to the localization length as well as interaction corrections to the 
inverse compressibility of the system [112] . Separating the two is not easy, 
while computationally it requires both a calculation of the sensitivity of the 
ground state to flux, as well as the dependence of the number of electrons in 
the system on the chemical potential. 
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Therefore, the study of the influence of interaction on the FO in the An- 
derson phase is not only interesting on its own account, but it establishes a 
new numerical method using a ground-state property which is convenient for 
a direct evaluation of the localization length for not too strong disorder. Us- 
ing this method we show that the localization length decreases as a function 
of the interaction strength, in correspondence to Eq. (14. ip . 

4.1.3 Chapter's Outline 

The rest of this chapter is organized as follows. In section 14.21 the model 
Hamiltonian and its diagonalization method are presented, followed by a 
description of the method used to extract the FO data from the ground- 
state wave functions. In section 14.31 an analytic calculation is preformed for 
the non-interacting case, resulting in an exact formula for the FO for a ID 
tight-binding model. Next we present the results for the two phases. The 
results for the TLL phase, which is replaced by an AI phase once disorder is 
introduced, are shown in section 14. 4[ The fact that the decay length can be 
used to approximate the localization length in the AI regime is also proven in 
this section. The results for the MI regime are shown in section H75| where we 
also provide an explanation to the difference found between the two phases. 
We conclude in section 14. 6[ in which we point out some of our planned future 
studies. 



4.2.1 Hamiltonian 

The system under investigation is composed of spinless electrons on a ID 
lattice coupled to an impurity in one end. We model the ID wire by a 
lattice of size L with repulsive nearest neighbor (NN) interactions and with 
an on-site disorder. The system's Hamiltonian is thus given by 
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(4.3) 
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where ej are the random on-site energies, taken from a uniform distribution 
in the range [—W/2, W/2], I is the NN interaction strength (/ > 0), and t is 
the hopping matrix element between NNs, henceforth taken as unity, cj (cj) 
is the creation (annihilation) operator of a spinless electron at site j in the 
lattice, and a positive background is included in the interaction term. 

Without the disorder term, a similar system - in the limit L ^ oo and 
with periodic boundary conditions - has a well known exact solution. De- 
pending on the interaction strength, the wire can be either metallic or insu- 
lating. As was detailed in chapter [H the metallic phase is described by TLL, 
occurring for / < 2t, and the insulating phase, in which I > 2t, is a CDW. In 
chapter [3] we have shown that wires of the order of a few hundreds sites lead 
to a similar phase diagram, even when employing open boundary conditions. 

Introducing an impurity at one end of the wire results in adding the 
following term to the Hamiltonian: 

Himp = eo4co - V{clci + H.c.) + /(cjco - ^)(clci - ^), (4.4) 

where eo describes the impurity strength, and V is the hopping matrix ele- 
ment between the impurity and the first lead site. Along this chapter we use 
eo > and V = t. 

The resulting Hamiltonian H = Hwire + -f^fmp describes a disordered ID 
wire of length L {1 < j < L), which is coupled to a single level at one of its 
edges (j = 0). Practically the j = site is equivalent to any other site, except 
for having a constant on-site energy, whereas the other sites have energies 
drawn from a distribution with a zero average over different realizations. 

4.2.2 Diagonalization Method 

The Hamiltonian H was diagonalized using the finite-size density-matrix 
renormalization group (DMRG) method [16], and the occupation of the lat- 
tice sites were calculated, for different values of eo, W and /. The size of 
the wire was up to L = 500 sites. During the renormalization process the 
number of particles in the system is not fixed, so that the results describe 
the experimentally realizable situation of a finite section of a ID wire which 
is coupled to a dot and to an external electron reservoir. Yet, the calcu- 
lated density remains close to half filling in all the calculated scenarios (even 
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in the presence of disorder) since the interaction term contains a positive 
background, and the calculation is done for = 0. 

4.2.3 Extracting the Friedel Oscillations Decay 

In the TLL phase, when no disorder is present {W = 0), Hwire has a particle- 
hole symmetry, and the particle density of the wire's ground state is fiat, 
with filling factor n = 1/2. In this case 2kp = n and the oscillating part 
of Eq. ( ]4.2p alternates according to (—1)-'. Denoting by (n™*"*^) 
the electron density at site j of the wire when coupled (not coupled) to the 
dot, one has nj^^'^ = n = 1/2 for any j. Clearly this is not the case in 
the presence of the impurity, and the effect of the impurity is measured by 
n. A typical result of Nj, showing the 2kF oscillations 



wire+imp 



caused by the impurity at j = 0, is presented in Fig. 14. 1[ 



Z 




Figure 4.1: A typical form of Friedel oscillations, without disorder and 
without interactions. The results shown are for eo = 10 and L = 280. The 
impurity is located at j = and the population of the first 20 lead sites is 
shown. 



When PF 7^ 0, on the other hand, although the average filling factor is 
still n ^ 1/2, there is no local symmetry between particles and holes, and 
the disorder effects are seen in the fluctuations of the electron density. The 
density oscillations generated by the additional impurity are then difficult to 
discern, since in a distance of a few lattice sites from the impurity the disorder 
fluctuations are dominant. A typical result of Nj together with N'j = nj^^'^—n 
(the electron density of the disordered wire without an impurity), is shown 
in the upper panel of Fig. 14.21 
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Figure 4.2: Typical FO for a disordered sample with L = 280, W = 0.1 
and eo = 10 (without interactions, TLL phase). The upper panel shows 
Nj (circles) and Nj (squares), and the lower panel presents the difference 
between them (ANj). The FO are observed much better using AA'^- instead 



ofiV,. 



However, the influence of the impurity can be observed by isolating the 
density fluctuations created by the disorder. This is achieved by comparing 
the electron density of the two cases shown in the upper panel of Fig. 14.21 i.e., 
one with the additional impurity and the other without it, for every disorder 
realization. Averaging over realizations is thus done for 

ANj = Nj - N^ = _ nJ*-, (4.5) 

instead of just averaging over Nj. The curve of ANj in the lower panel of 
Fig. 14.21 is for the same realization as in the upper panel. It is obvious that 
the FO which were hardly seen for Nj become clear once ANj is considered. 

In the CDW phase, when no disorder is considered [W = 0), the ground 
state of the system is twofold degenerate. This degeneracy is broken, however, 
once a pinning impurity, denoted by e^^^ — ^ 0+, is coupled to one end of the 
wire, and the particle density {Nj, for j = 1 . . . L) shows a, 2kp modulation 
(see chapter [3l). This is different from the TLL phase, in which the density 
without the dot is flat. 

When the pinning impurity is replaced by a dot level with eg ^ eQ^"*, the 
particle density in the wire (Nj) is changed by an additional oscillating 2k p 
term. One should notice that once the dot is coupled, a new CDW state 
emerges, having also 2kF oscillations, but with a different amplitude. The 
difference between these two states has a 2/c^ oscillation, with an exponential 
decay from its value at the edge of the wire. 
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Figure 4.3: Typical FO for a clean sample with L = 280 for a CDW 
with eo = 10 and / = 2.5. The upper panel shows Nj (circles) and Nj 
(squares), and the lower panel presents their difference ^Nj. 



In order to calculate the density difference between the cases when the 
quantum dot is coupled or uncoupled to the wire, it is easy to convince 
oneself that the definition of ANj in Eq. (14. 5 p should work also in the CDW 
case. A typical result of A^^- vs. Nj for a CDW state, and the resulting ANj, 
showing the 2k f oscillations caused by the dot orbital at j = 0, is presented 
in Fig. 01 

When W ^ 0, the CDW ground state is no longer degenerate, and the 
infinitesimal pinning impurity is not required. The disorder itself pins the 
CDW to different places on the lattice, with the ability to break the long 
range order of the clean CDW by localized solitons, with a density which 
depends on the disorder strength [91]. Yet, when a dot level with eo ^ is 
connected to one side of the wire, the local effect in its vicinity is stronger than 
the pinning caused by the disorder. This results in a change of the particle 
density near the dot, and this change decreases with distance. It turns out 
that the definition of ANj in Eq. (14.51) is suitable for the disordered case as 
well, since it cancels out the disorder pinning effects which are the same for 
the two cases, isolating the density fluctuations created by the dot. 

A typical picture of ANj for a disordered CDW sample is presented in 
Fig. 14.41 Whereas the upper panel shows the density of the two similar 
systems, one which is coupled to the quantum dot and the other is not, the 
lower panel presents the difference between these two densities, and the decay 
of the oscillations can be clearly seen. 
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Figure 4.4: Typical FO for a single disordered sample with L = 280, 
= 0.1 and eo = 10, for a CDW with / = 3. The upper panel shows Nj 
(circles) and (squares), and the lower panel presents their difference 
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4.3 Non-Interacting Case 

4.3.1 Exact Calculation of Friedel Oscillations 

We now wish to exactly calculate A^(m), the density of electrons in site m, 
of a half filled ID tight-binding lead, which is coupled to an impurity, in the 
asymptotic (m ^ 1) limit. The system is described by the Hamiltonian 

L-l 

H = eocjco - V{c{c^ + H.c.) - t ^(cjc,+i + H.c). (4.6) 

3=1 

N{m) can be calculated using the retarded Green function of an electron 
in the m'th site G^{uj] m, m) [70], and the relation (for a half filled band) 

1 r° 

N{m) = $5 / G^iuj;m,m)du, (4.7) 

where we are possibly neglecting bound states with energy lower than — 2t, 
which give exponentially small contributions for large m. The Green function 
itself is determined by 

G^{uj;m,m) = G^{uj;m,m) + (4.8) 
G^(cc;;m, 1) ■ V ■ G(u;;0,0) ■ V ■ G^{uj;l,m). 

In this expression G^(c(j;m, /) is the bare (i.e., without dot) lead Green 
function, while G{uj]0,0) = (tu — eo — T,{lj))~^ is the dot's Green function. 
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where ^(to) = ^{^t ~ — is the self energy of the dot, which was 

calculated in chapter [31 The first term in the RHS of the equation simply 
gives the constant n = 1/2 occupation in the absence of the dot. 

Substituting the known wave functions and energies of the tight-binding 
Hamiltonian one finds 



^D, a cos(A;a(m — /)) — cos(A;a(m + /)) 
G^{u]m,l) = - 2^ 



fc>0 



uo + 2t cos{ka) 



(4.9) 



where k = j^Uk, for integer n^. Transforming to integration over unit circle 
in the complex plane leads to 



-- + mA-(-)^ 

2t V 2t' 



(4.10) 



Combining Eqs. (gZD, dM]) and fl4.10p . one can get, 



AN{m) = N{m) - 1/2 = -^^^ 



duj- 



2m 



and by substituting oj = —2tcos9, we find 



AN{m) = — : 



w/2 



de- 



sin(6')e 



i2md 



T^ti J_^/2 2tcos(^) + Co - ^e^'^' 



(4.12) 



One now defines z = e *^ in order to get 



AN{m) 



2^ 



dz 



-2m 



z z'^ + e^z/t + l-V^/t'^'' 



(4.13) 



where the integration is over the right half of the unit circle, between the 
points ±1 on the imaginary axis (contour A in Fig. 14. 5p . 



<(C 



Figure 4.5: The integration contours A and B-C-D which connect the 
points [0,-1] and [0,1]. 
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Next we deform our contour to the contour B-C-D in Fig. 14.51 In doing 
so we neglect the contribution of poles which may occur inside the closed 
line A-B-C-D. These represent states bound at the impurity, and as we have 
mentioned above, contribute exponentially small terms for large m. The 
integration in parts B and D is done by defining z = ±ix, respectively, 
X G [1, oo), while the contribution of the semicircle C vanishes as its radius 
goes to infinity. Therefore we get 

For m 3> 1 the term x~2"^-i varies much faster than the other terms, and 
the rest of the integrand can be evaluated at x ~ 1 to give ^2/^2^^^^/^ - One 
thus gets the final form 

^N(,n)^t}r^M^Yly\ (4.15) 

Tim 1/ 2 

4.3.2 Numerical Results 

In order to check the accuracy of the FO data obtained by the DMRG 
method, we begin by a comparison of the numerical results of the clean 
non-interacting case, to the exact formula we've just calculated. For our 
model, substituting V = t'm. Eq. (14.151) . one gets the amplitude of the FO as 

A{eo) = --(t,/t+\\ . (4.16) 



7r 



As expected, the amplitude A does depend on eo, and this dependence 
is presented in Fig. 14.61 Curves obtained for different values of eo have the 
same slope, but not the same amplitude. The limits of eo ^ and eo ^ cxd 
are well understood, because in both of them the impurity does not play 
any role, the lead has a hard wall boundary, and the particle-hole symmetry 
imposes that the FO amplitude goes to zero. For finite values of eo, the 
behavior of the amplitude shown in the inset of Fig. 14.61 is compared to the 
exact relation Eq. (14.161) . The correspondence between the numerical results 
and this formula is excellent. 
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Figure 4.6: The FO decay (in log-log scale) for different values of eo 
(shown next to the curves) for clean samples without interaction. The 
slope - representing the decay exponent - is constant, and the only effect of 
eo is a change of the amplitude A. The curves for eo = 2.0, 10.0 are drawn 
with symbols. Inset: the dependence of A on eo for V = t together with 
the exact formula Eq. ()4.16p . 



4.4 From the TLL Phase Towards Anderson 
Insulator 

4.4.1 A Clean Lead: Extracting the TLL Parameter 

We now move to the interacting case. For < / < 2t, i.e., when the fermions 
in the lattice are described by the TLL theory, the decay is expected to be 
proportional |109] to In our model the TLL parameter k is given by [77] 

(4.17) 



2cos-i[-//(2t)]' 

For non-interacting particles one gets k = 1 so that the oscillations decay 
as while in the interacting regime a monotonic decrease of k toward 
the limit k, = 1/2 occurs as a function of interaction strength. Thus, as / 
becomes stronger, n decreases, and a slower decay is predicted. This trend 
is seen in the DMRG results presented in Fig. 14.71 

In the inset of Fig. 14.71 the results obtained for k, by fitting the FO 
decay of a 500 sites wire, to the predicted decay of x~'^, are presented to- 
gether with the theory prediction for k,{I) of Eq. fl4.17p . As can be seen, 
the results are in good accordance with the theory for interaction strength 
I /t < 1. Similar results were obtained using other implementation of the 
DMRG method (with a constant number of particles) |110] . and by func- 
tional renormalization-group studies |111] . In these works it was argued that 
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Figure 4.7: The DMRG results for the FO decay in log-log scale for 
/ = 0,0.5, 1 and 1.5 (bottom to top), which correspond to the TLL phase, 
with eo = 10 and L = 280 and without disorder. As / increases, the decay 
gets slower. Inset: The interaction parameter k as found by fitting the 
FO decay to x~'^ (symbols), together with the theory prediction Eq. ()4.17p 
(line). The DMRG results were obtained for L = 500 with eo = 1 (circles) 
and eo = 10 (stars). 



for the system sizes treated, the asymptotic regime in which the x~'^ behav- 
ior is predicted is not yet reached. In Ref. |lll j it was shown that using the 
functional renormalizat ion-group method, which is argued to be as accurate 
as the DMRG method, even L of the order of 10^ is not sufficient to obtain 
the values of k of Eq. KT7\ for J/t > 1. 

4.4.2 A Disordered Lead: Extracting the Localization 
Length 

We now turn on the disorder by taking W 0. In this case the results of AAT,- 
are averaged over 100 different realizations of disorder, which are sufficient to 
reduce the sampling error to less than one percent. In Fig. 14.81 the averaged 
particle density for W = 0.1 is shown and compared to the W = case 
for various interaction strengths. As can be seen, for small values of the 
interaction the effect of disorder is very weak, while for large values of /, the 
FO decay faster in the presence of disorder. Zooming into these curves, it 
can be shown that the effect of disorder is to multiply the clean FO decay 
by an exponential factor e~^^^, where ^ is a characteristic decay length. 

For each strength of the interaction, one can rescale the disordered W 
curves, to the clean W = one by simply multiplying it by e^^^, using ^ as a 
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Figure 4.8: The decay of FO for I = 0.0, 0.5, 1.0 and 1.5 (bottom to top) 
with L = 500 and eo = 10. The symbols are for W = 0.1 (AI phase) and 
the hnes are for the clean (W = 0) case (TLL). The disorder effect becomes 
significant for large values of / where the localization length is small. The 
average was done over 100 realizations. 

fitting parameter. As can be seen in Fig. 14.91 by using this rescaling method, 
the averaged disordered data collapses on the curves of the clean sample. 

The dependence of the decay length ^ on the interaction strength I is 
shown in the inset of Fig. 14. 9[ We shall now show that this quantity is 
effectively the mobility localization length. 

The effect of disorder in the continuum limit can be divided to the forward 
and backward scattering terms. Whereas the backward scattering term is 
related to the conductance and thus to the localization length of the electrons, 
forward scattering processes contribute only to the decay length of the FO, 
but not to localization. Thus, at first sight ^ does not necessarily correspond 
to the localization length. Nevertheless, in this case one can argue that the 
contribution of the forward scattering process to ^ is small and therefore ^ 
is a good measure of the localization length. 

Using standard bosonization technique it can be shown that the forward 
scattering processes result in the following term in the Hamiltonian: 



where is the TLL field which is related to the density operator by p{x) = 
— :^ V 0(2^) cind ri{x) is the g ^ component of the random potential. Since 
the TLL Hamiltonian (see chapter [H Eq. fll.30p ) 




(4.18) 




(4.19) 
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Figure 4.9: The rescaled decay of FO for the W = 0.1 curves over 
the clean curves of the TLL phase with different interaction strengths. 
Inset: the localization length found by the best fit for each value of / for 
W = 0.1 (circles) and W = 0.2 (squares). The lines correspond to the 
theory prediction of Eq. (|4.1|) with = 7000 and 1500 respectively. 

depends on (p only through (v0(a;))^, it is easy to show that by a redefinition 
of the field = — ^/^ dyr]{y) one can incorporate the Hfs term inside -^tll 
and get a similar form of Hamiltonian. Therefore, the forward scattering term 
is not expected to change the physics of the system. 

Nevertheless, it was shown that this redefinition of the field has an effect 
on the correlation functions [31J. This results in a decay of the density-density 
correlation function, which is, practically, the quantity we measure, and this 
decay is not related to the conductance. It is an exponential decay of the 
form e""^/', where / = 2^(k)^' -^f forward scattering strength of 

the disorder (defined in the non-interacting case). 

For the decay described by the characteristic length /, one can find, using 
the Bethe ansatz solution, the factor u/k, for each value of I. It is easy to 
show that u/k, and thus I, are monotonically increasing functions of J, as 
opposed to the FO decay length (see Fig. 14.91 in the inset). 

Moreover, one can estimate / quantitatively for the system we consider. 
The factor u/k, found from Bethe ansatz solution ranges from u/k = 2 for 
I = to u/k, ^ 4.5 for I = 1.5t. Denoting the amplitude of the disorder 
correlation function by D, i.e., {V{x)V{x')) = D6{x — x'), one finds that Df 
and (the forward and backward scattering disorder strengths) are of the 
same order of magnitude as D. For non- interacting spinless electrons in a 
ID lattice [m] 1/Db ^ 100/W^. Substituting W = 0.1, one gets / of the 
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order of 10^, which is much longer than the observed decay length. 

We thus conclude that the backward scattering processes are much more 
significant in the model treated, thus ^ is a very good approximation to the 
localization length, and its interaction dependence should be described by 
Eq. dH]). 

Using the prediction of Eq. (14.11) with the value of the disorder we employ 
along this chapter (order of 10~^), and recalling that without interactions 
^0 ~ 100/W^, the localization length should range between ^(/ = 0) ~ 
10^, which is much larger than the lattice sizes we considered, and thus 
almost doesn't influence the electron density, to .^(J = 2) ^ 10^, in which the 
disorder effect should indeed be much more dominant, in agreement with the 
qualitative results presented in Fig. 14.81 

The quantitative data shown in the inset of Fig. 14.91 flts the theoretical 
predictions of Eq. (14. ip for not too weak interactions. For weak interactions 
(/ < 0.5) no such flt was found, which is however expected, since for this 
regime the theoretical localization length is much larger than the wire length. 
The fact that the best flt to Eq. (jH]) was for ~ 7000 (1500) for W = 0.1 
(0.2), and not the expected ~ 10000 (2500), can be attributed to the 
same reason, as well as to the neglected forward scattering term which is 
stronger for weak I. We also note that the exact choice of the wire slices 
over which the flt is done, can change slightly the values of ^. This, however, 
does not change the qualitative results, showing a monotonic decrease of ^ 
as a function of the interaction strength. 

To summarize, the effect of disorder on the FO decay in the Anderson 
regime can be described by an extra exponential decay of the FO, which 
depends on the localization length, of the form 

(AiV,) = A{-iyr^eM-j/m). (4.20) 

where the localization length ^(k) decreases monotonically as the interactions 
increase. 
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4.5 From The CDW Regime Towards Mott 
Insulator 

4.5.1 A Clean Lead: Extracting the CDW Correlation 
Length 

Since the CDW is an insulating phase, the decay of the 2kF oscillations 
without disorder is supposed to be exponential and the characteristic length 
is the CDW correlation length [91], i.e., oc exp{—x/(). In Fig. 14.101 such an 
exponential decay of AAT,- is shown on a semi-log scale for various interaction 
strengths. An exact Bethe ansatz solution [32] of our model gives the relation 
between the correlation length and the interaction as 

C~exp(7r/v///(2t)-l). (4.21) 

The correlation lengths extracted from the DMRG results are presented 
with a fit to the exact formula in the inset of Fig. 14.101 As can be seen, for 
/ not very close to the TLL-CDW transition point (which occurs at / = 2t), 
the results fit the theory very well. 




Figure 4.10: The oscillations decay in the CDW regime for various 
interaction strengths and without disorder (note the semi-log scale). As 
the interaction increases, the correlation length decreases and the decay is 
faster. Inset: the inverse correlation length of the CDW state for various 
interaction strengths (symbols) fitted to the theory prediction Eq. (j4.2ip . 

4.5.2 Decay Length in the Disordered CDW Phase 

For W ^ 0, ANj is averaged over 100 realizations, for which we expect a 
sampling error of the order of one percent. Assuming that the disorder adds 
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another exponential term to the oscillations decay, which is thus proportional 
to exp{~x / ( — X / C,) , there are two competing length scales - the decay length 
due to disorder (^) vs. the correlation length {(). For strong interactions 
and weak disorder C <^ ^ so that the disorder effect is hardly seen, but 
increasing the disorder or decreasing the interaction strength should result in 
a combination of the two exponential decays. The DMRG results, presented 
in Fig. 14. m show the disorder effect on the oscillations decay. For / = 2.5 
and / = 3 one can see faster decay for the disordered samples with W = 0.1. 
For stronger interaction larger disorder is required in order to affect the decay. 




10 20 ; 30 



Figure 4.11: The decay of the oscillations of a disordered CDW with 
I = 2.5,3 and 3.5 (top to bottom, note the semi-log scale). The lines 
correspond to the clean sample result, and the symbols to the averaged 
disordered data. For 1^ = 0.1 (circles) the disorder effect is clearly seen 
for / = 2.5 and / = 3 but not for / = 3.5 in which ^ is much larger than 
the correlation length For W = 0.2 (squares) ^ is small enough to affect 
the decay even for / = 3.5. Inset: multiplying AA'^- by e^^^ collapses the 
disordered data on the clean curves. 

Similarly to the AI phase, the extra decay length can be extracted by 
fitting, for each value of J, the W ^ curve multiplied by e^^^ to the W = 
one. Such a rescaling is presented in the inset of Fig. 14. Ill 

4.5.3 The Decay Length Dependence on Interaction: 
AI vs. MI 

As can be seen in Fig. 14.121 the decay length extracted for the disordered MI 
regime increases as a function of the interaction strength (for 2t < / < 3.5t), 
an opposite behavior to the AI case (/ < 2t) . Results for stronger values of I 
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Figure 4.12: The decay length due to disorder (^) in the TLL (/ < 2t) 

and in the CDW (/ > 2t) phases as a function of the interaction strength. 
Inset: zoom into the CDW regime. 

are not shown, since for too strong interactions the correlation length is very 
small, and thus the estimate of ^ is less accurate. 

These results reveal that as the interaction strength increases in the MI 
phase, the disorder effect decreases. In the AI phase, on the other hand, 
the disorder effect is enhanced with increasing interactions. The difference 
between these two behaviors results from the difference in the ground states of 
the two phases in the clean case. In our model there is a competition between 
the kinetic energy (the hopping term) and the potential (the interaction). 
The hopping term prefers the existence of a flat particle distribution whereas 
the interaction term prefers a CDW-like form. For different values of / the 
results of that competition are different: for I < 2t (the TLL phase) the 
hopping term wins, and the distribution is flat, while for I > 2t (the CDW 
phase) a CDW starts to form. 

Inside the clean TLL phase, as / increases, the CDW fluctuations are 
stronger. Yet, the average density proflle in the ground state remains flat 
because of the hopping term. But when disorder is introduced, the flat 
density state becomes less favorable than a state with a fluctuating density, 
the latter being preferred by both the disorder and the interactions. For a 
constant disorder, as the interactions become stronger, these fluctuations are 
enhanced, so the disorder effect increases. 

In the CDW phase, on the other hand, without disorder, the interaction 
wins over the hopping, and the ground state has a CDW form. Turning on the 
disorder might change the particle distribution, e.g. by allowing an electron 
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to move into a site with lower on-site energy, but this results in raising the 
interaction energy. As the interaction strength gets stronger, the probability 
of such a process decreases, so that the actual effect of the disorder is getting 
weaker. 

4.6 Conclusions and Future Prospects 

In conclusion, we have shown that the FO envelope is affected by both the 
interaction strength, and the disorder strength. In the TLL phase and the 
resulting disordered AI phase, interactions actually enhance the effect of FO 
since it drops with a weaker power law while disorder decreases the 
FO since it adds an exponential factor to its decay form. We have shown 
that the length scale for this exponential decay is a good approximation to 
the mobility localization length, since it is only weakly influenced by forward 
scattering processes for weak disorder. Thus we established a convenient way 
to evaluate the dependence of the localization length on disorder and inter- 
action using only the ground-state properties of the system. Qualitatively, 
the localization length as a function of interaction for a given weak disorder 
always decreases. As long as the localization length is not much longer than 
the wire's length, the localization length behavior is quantitatively described 
by the renormalization-group results [88] . 

However, while the decay length of the 2k p oscillations envelope due to 
disorder is monotonically decreasing in the AI phase, we have shown that it is 
monotonically increasing in the disordered MI phase. The difference between 
these two regimes is explained by the difference between the ground states of 
the clean samples in each case. In the AI phase the pure ground state is flat, 
and both the disorder and the interactions try to introduce fluctuations in 
it. In the MI phase, on the other hand, the pure ground state oscillates with 
a 2kF wave vector, and these oscillations are enhanced by the interactions 
and reduced by the disorder. As a result, the disorder effect (for a constant 
disorder strength) is getting weaker, in the MI phase, as the interactions are 
enhanced. 

In addition, we have analytically described the dependence of the FO 
amplitude on the impurity strength and on the dot-lead coupling, in the 
non-interacting case. 
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* * * 

During the discussion of the disordered MI phase, we have mentioned 
that the strength of the disorder has a significant impact on the phase of the 
wire. While for a weak enough disorder the wire can still be described by a 
MI phase, for a strong disorder it turns to an AI phase. A few recent studies 
[nH ES] point out the possibility that another phase exists between the MI 
and the AI, which is the Mott glass (MG) phase. 

The MG phase is characterized by partial features of both the AI and the 
MI phases. For example, it does not have a long range order, since it includes 
ordered regions which are separated by some kind of domain walls (solitons 
in the CDW language). In addition, the MG energy spectrum is gapless. In 
both these properties the MG phase is thus similar to the AI one. On the 
other hand, the low-lying excitations are not charged, but describe a motion 
of the domain walls, and as a result the MG phase is incompressible, like the 
MI phase. 

The model used in Refs. [Ml ES] consists of spinless electrons in a disor- 
dered potential, as in our case. Nevertheless, in order to assure the existence 
of a clean MI phase, as opposed to the TLL phase, it introduces a2kp com- 
mensurate potential as well. Therefore, when there is no disorder, a MI 
phase is obtained for the entire range of interactions. The effect of disorder 
is now added to that of the periodic potential. The resulting phase diagram 
shows that for a non-vanishing disorder strength, a gradual increase of the 
commensurate potential causes the wire to turn from an AI (when there is 
no periodic potential), via a MG, towards a MI (strong periodic potential). 

It is therefore clear that a finite regime of the MG phase can be obtained 
using a similar model to the one we use, by adding a periodic potential to 
the Hamiltonian. An interesting question can then be raised: How will the 
decay of the FO look like? Can the shape of the FO decay help to identify the 
borders of the MG phase? The results of this chapter show that one is able to 
differentiate between the AI and the MI phases by looking at the dependence 
of the decay length on the interactions. In the case of the MG phase, such 
an identification is much less clear. However, the versatile properties of the 
MG phase might be helpful for that task, which we leave for a future study. 
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Finally we remark to the experimental relevance of this work. The the- 
oretical treatment of disorder usually involves statistics over an ensemble of 
many samples which is usually hard to obtain experimentally. Furthermore, 
in the case we deal here, a measurement of FO on a disordered sample seems 
at first sight daunting. However, the simple method we suggest in order 
to deal with the disorder, is in principle experimentally feasible, and solves 
these two difficulties. 

Once a technical method for measuring the electron density is established, 
it should be used twice for each sample, before and after the coupling of the 
wire to the dot. In principle, by using a gate it should be possible to eliminate 
the coupling between the dot and the wire. Our results, as can be seen in 
Figs. 14.21 and 14.41 which present typical results for a particular realization, 
point out that the difference between these two measurements should show 
a very clear FO, even for a specific sample. 



Chapter 5 

A Disordered QD with 
Interacting Electrons 

In this chapter we explore some properties of a two-dimensional quantum dot, 
without considering its coupling to external leads. We model the quantum 
dot by a two-dimensional lattice with either nearest neighbor or Coulombic 
interactions, which is occupied by spinless particles. Since the size of the 
relevant Hilbert space is huge, such a problem cannot be solved by an exact 
diagonalization (except for very small systems), and we thus use the particle- 
hole density-matrix renormalizat ion-group (PH-DMRG) method, which was 
detailed in section 12. 5[ 

One of the most important quantities involved in any physical system is 
the energy spectrum, and especially the ground-state energy. In this chapter 
we show that the PH-DMRG method can be used in order to approximate 
the ground-state energy of the disordered quantum dot with interacting elec- 
trons, and that it leads to results which are much more accurate than those of 
the Hartree-Fock approximation. Moreover, following a comparison to other 
approximate methods, we suggest an improvement of the PH-DMRG trun- 
cation algorithm, which reduces the error rate of the traditional^ PH-DMRG 
by almost 30 percents. 

As an application for the improved PH-DMRG method we calculate the 
ground-state energies for different numbers of electrons, and find the addition 
spectrum of the system. We compare the PH-DMRG results to those of the 

^ It's a bit weird to call it traditional since the known number of PH-DMRG applications 
can be count currently using one hand. However, by "traditional" we refer to the method 
as it is documented in [60j and |63j . 
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Hartree-Fock approximation in three aspects: the error rate, the average and 
the fluctuations of the addition spectrum. 

5.1 Introduction 

In the last decade there is a growing interest in the low-temperature physics 
of disordered many-particle systems, such as electron dephasing due to in- 
teractions [114j and two-dimensional (2D) 'metal- insulator' transition |115j . 
Transport properties through quantum dots (QDs) have also been recently 
investigated and shown to exhibit some interesting phenomena in the pres- 
ence of both interactions and disorder [3]. An analytical treatment of these 
problems is difficult, unfortunately, since both the disorder and the interac- 
tions cannot be considered as a small perturbation. A traditional numerical 
treatment of such problems is restricted to small systems, since the size of 
the many-particle Hilbert space grows exponentially with the system size. 

During recent years, a few methods were established in order to decrease 
the size of the Hilbert space to a size which is computationally feasible. One 
way is to define an iterative order in which the system is treated, and use a 
sophisticated truncation method between the iterations to reduce the space 
size. This is the idea behind the ensemble of renormalization group methods, 
among which the density-matrix renormalization-group (DMRG) method is 
an honored member (see section [27il) . 

A different approach uses a predefined constraint which is checked before 
states are inserted into the Hamiltonian. The entire system may be treated 
immediately, yet not all the system states are taken into the Hamiltonian. 
Therefore, the matrix size one needs to diagonalize is smaller than the entire 
Hilbert space dimension, and hopefully small enough to be exactly solved. 
Yet, if the constraint is defined on the many-particle space, then although 
the final matrix size might be small, the calculation time required grows 
exponentially with the system size. Thus such a constraint does not solve 
the main issue. 

However, Berkovits in Ref. |116] has investigated constraints which are 
checked against the single-particle states before they are used to build the 
multi-particle basis. The usage of this approach was demonstrated on a 
ground-state energy calculation of a disordered 2D lattice of 4 rows and 6 
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columns containing 10 interacting electrons, for which the full Hilbert space 
contains 1, 961, 256 states. The error rate, or the discrepancy, can be defined 
a.s D{x) = {\x' — x\ / \x\) , where x is an exact quantity, x' is an approximation 
for X, and the average is done over different realizations of the disorder. By 
calculating the value of D{Egs), one thus has a good estimate about the 
accuracy of the approximation method. 

In this chapter we quote the results obtained in Ref. |116j by two meth- 
ods, which we denote as "energy-cut" and "generation-cut". In both methods 
the first step is to execute a self-consistent Hartree-Fock (HF) approxima- 
tion, and obtain a sorted "list" of wave-functions and energy levels, which 
partly incorporate the interaction effects. In the first method, the energy- 
cut, one neglects the states with the highest single-particle energies. Such 
states should, intuitively, have the smallest contribution to the many-particle 
ground-state energy. For a system with sites and Ue electrons, the full 
Hamiltonian matrix is of size M = (^J . Using the energy-cut method, one 
takes only the lowest m/j single-particle states, and thus diagonalizes a ma- 
trix of size Mr = (^^) ■ In the results of Ref. |116j , the energy-cut method 
was used for the 4x6 lattice with 10 electrons, by taking up to = 18 (out 
of the 24) HF states, which leads to a maximal matrix size of M/j = 43, 758. 
The best discrepancy obtained by this method was D{Egs) ^2.5 percent, a 
significant improvement of the HF results, for which D{Egs) ^ 4.5 percent. 

The second method, generation-cut, has obtained better results than the 
energy-cut method. This method is based on the localization of the Fock 
space. Since the interaction term is a two-body operator, only many-body 
states which differ by at most two electron-hole pairs are coupled by the 
Hamiltonian. It was proven pj.7] that the average contribution of a state 
containing k electron-hole pairs to the exact ground state is proportional 
to exp{—k/^p), where C,f is the Fock space localization length. Considering 
also the number of states in the k-th electron-hole generation, (^^""') {"'j^) , 
one finds out |116] that the weight of generations falls off exponentially as 
long as > \n[{N - ne)nj{k + if]. 

Therefore one can consider in the approximated Hamiltonian only states 
with a small number of particle-hole pairs. For exactly / particle-hole pairs, 
the number of states is j""") {y) ■ In the results of Ref. |116] , the generation- 
cut method was used by taking up to 3 particle-hole pairs, which leads to a 
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maximal matrix size of 47,916. The lowest error obtained by this method 
was D{Egs) ^ 1.5 percent. 

Unfortunately, one cannot simply 'cut and paste' these couple of methods 
for larger systems. As the size of the system increases, the number of impor- 
tant HF states (for the energy-cut method), and the number of important 
particle-hole generations (for the generation-cut) are expected to increase. 
Moreover, even for the same number of chosen HF states and particle-hole 
generations, the matrix size is exponential in the system size, and it becomes 
soon too large to be solved exactly. Therefore, a treatment of larger sys- 
tems will experience similar problems to those of the exact diagonalization 
technique. 

Thus there is a great need for an alternative method, which can give 
an accurate approximation to the ground-state energy of a 2D QD, yet can 
be extended to larger systems. In this chapter we check if the PH-DMRG 
method is suitable for that task. 

5.1.1 Chapter's Outline 

The rest of the chapter is organized as follows. In the following section 
we briefly describe the Hamiltonian model and the numerical PH-DMRG 
method, which was extensively discussed in section 12.51 In section 15.31 our 
main results are presented. We check the accuracy of the PH-DMRG method 
in approximating the ground-state energy of an ensemble of disordered lat- 
tices with interaction. We show that the PH-DMRG method can get an 
accuracy which is similar to that of the energy-cut technique, whereas the 
generation- cut method obtains better results. Analyzing some disadvantages 
of the PH-DMRG algorithm, we then suggest an improvement of the PH- 
DMRG truncation method. This leads to an improvement of more than 30 
percents of the PH-DMRG accuracy, making it comparable to the generation- 
cut technique as well. 

In section 15.41 we compare the PH-DMRG accuracy between different 
schemes: Short-range vs. long-range interactions, and intermediate vs. strong 
interactions. 

As an additional application we show in section 15.51 the results for a 
calculation of the addition spectrum using PH-DMRG, when either nearest 
neighbor (NN) or Coulombic interactions are present. We show that in both 
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cases the accuracy of the PH-DMRG method is much better than that of the 
HF approximation. 

We conclude the chapter in section \5M in which we discuss the possibihty 
to use the PH-DMRG for larger systems, and we point out an optional future 
research. 

5.2 Model 

The system we treat, as a QD model, is composed of Ue spinless electrons on 
a 2D disordered lattice of A columns and B rows. The electrons can hop from 
one lattice site to one of its NNs, and either NN or Coulombic interactions are 
considered. Therefore, the Hamiltonian describing the physics of the system 
can be written as H = Hq + Hint, where Hq is given by Eq. fl2.14p . and Hint, 
for NN interactions, is given by Eq. (12.151) : 

^0 = ^ e^aLom - ^ (am^n + H.c), (5.1) 

m (jn,n) 

Hint = ^ oLalanClrn- 

{m,n) 

Here, a]^ and denote creation and annihilation operators of an electron 
in lattice site m, where (m, n) represents NN sites m and n. is the on-site 
energy of site m, chosen randomly from a uniform distribution [—W/2, W/2], 
t is the hopping matrix element between NNs, and V is the NN interaction 
strength, t is conventionally taken as the energy unit. 

The ground-state energy of the system is calculated using the PH-DMRG 
method, which was described in details in section 12. 5[ The general idea is 
to divide the energy levels to those above (particle-states) and below (hole- 
states) Ep, and then treat these states iteratively. Starting from the vacuum 
state, in which all levels below Ep are filled, and all the others are empty, we 
add in each iteration one particle-state and one hole-state, starting from Ep 
and proceeding in both directions. We diagonalize the superblock composed 
of the states we already added, maintaining the number of particles constant. 
Each iteration ends with a truncation of the Hilbert space of both the particle 
block and the hole block, using their corresponding density matrices, after 
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which a new couple of states can be added again. The iteration process stops 
after all of the states were addec|^. 

As was emphasized in section l275l the accuracy of the PH-DMRG method 
depends mostly on the number of states, p, that are kept between successive 
iterations. In other words, in each iteration, p eigenvectors of the density 
matrix are taken, while all the others are neglected. Except the accuracy, 
p influences dramatically the computational resources required during the 
process. 

In general, as the Hamiltonian is more complicated, and as the number 
of both particle-states and hole-states increases, the size of the superblock 
Hamiltonian, whose creation and diagonalization are the main bottle-neck of 
the algorithm, is enlarged. Therefore, for different Hamiltonians, and even 
for different number of particles in a given model, the upper limit of p can 
change. In this study the maximal number of block states that we keep 
is p = 120, which is of the order of the number of states used in typical 
PH-DMRG applicationf§. 



5.3 Ground-State Energy Calculation 

We begin by presenting the PH-DMRG results for the calculation of the 
ground-state energy in an interacting disordered dot. The Hamiltonian for 
such a model is given by Eq. (15. ip . We investigate the case of 10 interact- 
ing electrons in a disordered 4x6 lattice, with the same ensemble used in 
Ref . |116j . In this section we restrict ourselves to the case of NN interactions, 
with strength V = 3t, and we set the disorder strength to he W = 5t. 

Typical results, for a specific realization, are shown in Fig. 15.1( a). as a 
function of p, the number of block states kept. The results are compared to 
the HP results (green line), to the results of the energy-cut method, in which 
the lowest 18 HP states were used (magenta), and to those of the generation- 
cut method, in which up to 3 particle- hole generations were kept (brown). 
The exact results are also drawn (a red line). As p increases, the PH-DMRG 

comprehensive description is presented in section [2.51 
^Pmax ~ 100 was used in most implementations. However, for simpler models, which 
lead systematically to small superblock Hamiltonians, PH-DMRG was utilized withpmax ~ 
400. It should also be noted that the presence of disorder in our model requires also 
averaging over realizations. 
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approximation is better. Typical values of p in which the PH-DMRG method 
is already more accurate than the energy-cut method are 40 — 50 (for the 
sample shown p > 30 was sufficient). In rare samples larger values of p were 
required in order to improve the energy-cut results. However, a very slow 
convergence is seen, and the generation-cut method obtains the best results 
for up to p = 80, for all samples. 

Averaging over realizations makes the picture more clear, as one can 
see in Fig. 15.1( b). The average accuracy of the PH-DMRG calculation is 
improved very slowly by increasing p. The best PH-DMRG results shown 
(keeping 80 states) gets lower discrepancy than the energy-cut method, yet 
the generation- cut method obtains much more accurate results. 
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Figure 5.1: PH-DMRG results for the ground-state energy calculation 
as a function of p, the number of states kept during the renormalization 
process. The results are compared to various approximated results (see 
text), (a) The ground-state energy of a specific realization. The exact 
result is presented by the red line, (b) The discrepancy D{Egs) = {\E'gg — 
Egs\/\Egs\) averaged over an ensemble of 100 realizations. 



The comparison between the PH-DMRG results and those of the generation- 
cut method raise the following question: Can one find a way to improve the 
PH-DMRG process (besides the option to increase the block size)? In the 
following we present such a possible improvement, motivated by an analysis 
of the PH-DMRG truncation method. Recalling the DMRG principles, the 
main difference between the DMRG and previous numerical renormalization 
group (NRG) methods, is the truncation algorithm. In the NRG method, 
the truncation of states is based on their energies, while in the DMRG it 
is based on the density-matrix eigenvalues. The density-matrix eigenvectors 
with the highest eigenvalues are considered as the most important, and the 
rest of the states are neglected. 
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Nevertheless, the difference between these criteria, i.e., the NRG's lowest 
energy and the DMRG's highest density-matrix eigenvalues, is not the main 
reason for the DMRG success. This success originates mainly from the phys- 
ical situation to which the density matrix emulates. While the NRG trunca- 
tion method is based only on the sites which were already iteratively added, 
in the DMRG algorithm other sites are also included in the superblock. The 
superblock is composed of a "system" coupled to an "environment", which 
represents all the sites which were not yet included in the iteration process, 
and thus the truncation is based on a future prospect, which leads to better 
results. 

On the other hand, the superblock in the PH-DMRG does not consider 
any "future" states. The PH-DMRG process couples only the already-used 
particle-states and hole-states, but other states, those which were not yet 
inserted, are not part of the superblock. Therefore, the truncation does not 
take them into account, and this might limit its success. 

In order to give the truncation decision a wider point of view, one may 
use an additional condition to help choosing the states to continue with. 
Such a condition should thus take into account some further considerations. 
Based on the great success of the generation-cut in Ref. |116j . we try here to 
consider the number of particle-hole generations as an extra condition taken 
during the truncation step. 

As we have discussed above, the idea of the Fock space localization indi- 
cates that the weight of successive particle-hole generations decreases expo- 
nentially. Moreover, in a comparison between specific states at a given time 
during the process, those with smaller number of particle-hole generations 
are apriori more important. 

There are a few ways in which one can incorporate this criterion into 
the density-matrix truncation process. For example, when the states are 
weighted by their density-matrix eigenvalues, one can add a multiplicative 
factor for each eigenstate, based on the number of particle-hole pairs it con- 
sists of. Such a method should be accompanied by a systematic investigation 
of the weighting procedure, e.g. by a calculation of the relative probability 
for each number of particle-hole pairs according to the Fock space localization 
length, and we leave this subject for future studies. 
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Here we use a simpler condition. We define a maximal number k of 
preferred particle-hole generations, and states with f < k are kept in the 
first stage. If there are too many such states, we use the density-matrix 
eigenvalues, and take the states with the highest eigenvalues. After this first 
round, if there is still a room for more states, we take also states with f > k, 
according to their density-matrix order, until the maximal number of states 
p is reached. 

For the 4x6 system with 10 electrons, we've executed the suggested 
truncation method for k = 0,1, 2, 3, 4, and for each realization we have picked 
the lowest energy. However, it turns out that in 95 percents of the samples 
the k = 2 case resulted in the most accurate result, so that using a constant 
. ^ 2 would lead to s„„,.ar .esultfi F,g. O presents the same curves of 
Fig. 15.11 together with the results of the improved algorithm. In Fig. 15.2( a) 
the results of = 2 are compared to the traditional results for a specific 
realization, while in Fig. 15.2( b) the averaged discrepancy is shown. As can be 
clearly seen, the results obtained by the improved method, which are better 
by almost 30 percents than the regular PH-DMRG results, are comparable 
to the results of the generation-cut method of Ref. [116]. 
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Figure 5.2: The same as Fig. 15.11 together with the results of the im- 
proved PH-DMRG method (traditional version - black circles, improved 
version - blue crosses). Inset: zoom into the region of p > 60. Both panels 
show that the improved PH-DMRG method obtains similar results as the 
generation-cut method for p ~ 110 — 120. 



From these results it is quite clear that the PH-DMRG can, in principle, 
obtain the ground-state energy with an accuracy of the order of 1.5 per- 
cent, when the number of states p is of the order of 100. In the last section 
we discuss how these results are affected when larger dots are considered. 



^In the following sections we thus use k — 2. 
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We will conclude that the enlargement of the system has much more dra- 
matic influence for the generation-cut method, compared to the PH-DMRG. 
Therefore, the PH-DMRG is a good candidate for the task of approximating 
ground-state energies of disordered QDs. 

Nevertheless, we still need to check how sensitive the accuracy of the PH- 
DMRG method is when the interactions have long range, or when they are 
much stronger than those considered in this section. This will be discussed 
in the following section. 

5.4 Long- Range Interactions; Strong Interac- 
tions 

The Hamiltonian used in the previous section, as well as that of section 12.51 
utilized only short-range interactions. In the framework of the real-space 
DMRG method there is a huge impact when the interaction range increases. 
The DMRG iterations add subsequent sites one after another, and if long- 
range interactions are considered, much more data should be stored from 
previous steps. Practically, therefore, real-space DMRG applications tradi- 
tionally consider only nearest neighbors, or next nearest neighbors interac- 
tions. 

On the other hand, the interactions range does not affect the PH-DMRG 
method almost at all. It does of course change the first stages of the method, 
i.e., the single-particle diagonalization and the self consistent HF stage, be- 
cause the Hamiltonian is changed. For example, for Coulombic interactions 
the Hamiltonian is (compare to Hint in Eq- (15. ip ) 



where jr^ — r„| is the distance between sites m and n, measured in lattice 
units. Accordingly, the definition of the antisymmetric interaction element 
Va/s^s in Eq- (I2.20p should be modified to 




(5.2) 



V 




r. 



m 



- r. 



n 



{ <P*a{m)(p*p{n) [(p^{n)(p^{ni) - (ps{m)(p^{n 



)] 




0;(n)(/)^(m) [(t)^{n)(t)s{m) - </)5(m)(/)^(ra)] } 
I [4>*a{m)4>*p{n) - 0;(n)0^(m)] [0^(n)05(m) - (p&{m)(p. 
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However, except for that slight change, the entire iterations process, which is 
the most taxing stage of the PH-DMRG algorithm, remains the same. The 
bottom line is that since the iteration is done over states in momentum-space, 
it does not matter what the real-space range of interactions is. 

Therefore, the range of interactions does not change the feasibility of the 
PH-DMRG method at all. However, one should ask what effect does the 
interactions range have on the PH-DMRG accuracy. A related question is 
how the accuracy is changed when the interactions become stronger, and does 
the enhancement of interactions in the short-range and long-range cases have 
the same influence. 

In Fig. 15.31 we present the results of the discrepancy obtained for short- 
range (upper panel) and long-range (lower panel) interactions, as a function 
of the number of states kept, p. In both cases we compare the results for 
intermediate strength of interactions (black) and a strong one (red). These 
results were obtained by the PH-DMRG method for 4 x 6 lattice with 6 
electrons. Since the number of electrons is small (compared to the previous 
section), there are less hole-states, and the PH-DMRG accuracy is better. It 
should be noted that the same is true for the result of the HF approximation 
(shown in dashed lines). 
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Figure 5.3: The discrepancy D{Egs) obtained by the PH-DMRG calcu- 
lation as a function of the number of states kept, for an ensemble of 4 x 6 
disordered lattices occupied by 6 electrons. The accuracy is compared be- 
tween short-range (upper panel) and long-range (lower panel) interactions, 
and between intermediate interactions strength iy = Vc = 3t, black sym- 
bols) and a strong one iy = Vc = lOt, red symbols). The results of the HF 
are presented in dashed lines. Note the semi-log scale. 



For both short-range and long-range interactions, we find that the accu- 
racy is reduced when the interaction strength is enhanced. Therefore, as the 
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interaction strength increases, the number of states p which is required in 
order to obtain the same accuracy, increases. 

A more striking feature is the difference between the short-range and the 
long-range interactions in the discrepancy. For small values of p, the better 
results of the long-range case are explained by the fact that the HF approx- 
imation is known to work better for long-range interactions. Nevertheless, 
the improvement of the long-range results by the PH-DMRG method is fas- 
cinating, and the obtained averaged discrepancy for Vc = 3t with p = 90 
states is D{Egs) ~ 10~^, more than an order of magnitude better than the 
accuracy in the short-range case. 

We thus close this section with the following conclusions: The PH-DMRG 
can be used as an accurate method for calculating ground-state energies of 
disordered QDs with interactions. The physical parameters of the system 
can affect the accuracy significantly. The best accuracy is obtained for long- 
range interactions which are not too strong, when the number of electrons is 
not too large. 

5.5 Addition Spectrum Calculation 

As we have shown in the previous sections, the PH-DMRG method can be 
used in order to get accurate results for the ground-state energy of interacting 
2D systems. As a useful application we present in this section a calculation 
of the addition spectrum of a QD, accompanied by a comparison of the PH- 
DMRG results to those of exact diagonalization and of the HF method. 
The addition spectrum can be defined by (see section I1.1.3P 

A2 = Eg,{N,) - 2EgsiN, - 1) + Egs{N, - 2). (5.4) 

Therefore, for a calculation of the addition spectrum one needs the ground- 
state energies of 3 successive electron numbers for each realization. The 
results shown in the current section are for the ensemble of 4 x 6 samples 
used in the previous sections, occupied by 4, 5 and 6 electrons, with either NN 
or Coulomb interactions. In general, the results were better in the Coulombic 
case, because of a higher accuracy for each energy calculation. 

In Fig. I5.4l the results for the averaged discrepancy, D{A2), are shown, as 
a function of the number of block-states kept in the PH-DMRG calculation. 
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Figure 5.4: The averaged discrepancy D(A2) obtained by the improved 
PH-DMRG calculation, for a 4 x 6 lattice occupied by 4, 5 and 6 electrons, 
with NN (left panel) or Coulombic (right panel) interactions. The PH- 
DMRG results are shown as a function of p (symbols), together with the 
HF results (lines). Note the semi-log scale in the right panel. 



for both NN (left panel) and Coulombic (right) interactions. The HF ap- 
proximation obtains, in the Coulombic case, D{A2) ~ 1.7 percents. As can 
be seen, the PH-DMRG method, even with a very small number of states 
{p > 30), reduces significantly the error rate to a level of ~ 0.13 percents, an 
improvement of more than an order of magnitude. 

However, in the NN case the results are quite poor. The starting point of 
the PH-DMRG algorithm, i.e., the HF results, give an average error of almost 
20 (!) percents. The PH-DMRG improves it by a factor of 4, to the order 
of 5 percent, which is still a very high error rate. It is thus interesting to 
check whether the calculation of (A2) and 6A2 can give more accurate results. 
Typical results are shown in Fig. 15.51 ((A9) (left panel) and 6A2 (right panel), 
for NN interactions), and in Fig. 15.61 (for Coulomb interactions). In all cases 
the corresponding PH-DMRG results are compared to the exact solution 
and to the HF results. In the insets we present KAg) — (A2)|/|(A2)| and 
\5A'2-6A2\/\6A2\. 

As can be seen, in the NN case the PH-DMRG results for (A2) are very 
accurate for p > 30. On the other hand, p = 60 is still not sufficient in 
order to get accurate results of 6A2, which shows very slow convergence. In 
both cases, however, the PH-DMRG results are significantly more close to 
the exact results than the results obtained by the HF method. 

In the Coulombic case, the results for (A2) are similar to those of the NN 
case, and very small value of p is sufficient to get very accurate results. Notice 
that the results continue to fluctuate around the exact result. This results 
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Figure 5.5: The results of the PH-DMRG calculation for (a) (A2) and 
(b) 5A2 of the 4x6 system occupied by 4, 5 and 6 electrons with NN 
interactions, as functions of p, the number of states kept. The green lines 
correspond to the HF results, and the red lines to the exact ones. Insets: 
KA'a) - (A2)|/KA2)| and I^A'g - <5A2|/|5A2| as functions of p. 



from the fact that the approximation for A2 is done using 3 different ground- 
state approximations, which have different convergence rates. However, for 
all values of p, the error is less than 0.1 percent. For 6A2, on the other 
hand, the convergence with increasing p is clearly seen. As can be seen, the 
Coulombic case leads to much more accurate results than those of the NN 
interactions, giving an error rate of only 0.8 percents. 
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Figure 5.6: The same as Fig l5.5| but with Coulombic interactions. Note 
the semi-log scale in the insets. 




To conclude this section, we define the improvement-factor of the PH- 
DMRG by the ratio between the HF discrepancy, and that of the PH-DMRG. 
Table 15.11 summarizes the error rates obtained by the PH-DMRG (using 
p = 60) and the respective improvement factors (in parentheses) for the 
different cases we've treated. 

It is easy to see that for Coulomb interactions the PH-DMRG improves 
all results related to the addition spectrum by an order of magnitude, and 
leads to error rates of less than 1 percent for all cases. However, for the NN 
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NN 


Coulomb 


(A2) 
SA2 


0.054314 (3.6) 
0.000965 (55.1) 
0.055554 (3.4) 


0.001330 (12.6) 
0.000352 (16.0) 
0.007794 (9.2) 



Table 5.1: PH-DMRG error rates for A2 calculation 



interactions, a small error rate and a significant improvement factor are seen 
only for (A2), while modest factors are obtained for D{A2) and 6A2, with 
error rates larger than 5 percents. 

5.6 Conclusions and Future Prospects 

In conclusion, we have seen in this chapter that the PH-DMRG method can 
be used for a calculation of the ground-state energy in disordered systems 
with interactions. We have analyzed some disadvantages of the method and 
suggested an improvement of the traditional implementation. We have also 
compared the accuracy of the ground-state calculation between long-range 
and short-range interactions, and between intermediate and strong interac- 
tions. We have found that the PH-DMRG works better when the HF approx- 
imation is better. The best accuracy is obtained for long-range interactions 
which are weak or intermediate, when the number of electrons is not near 
half filling. 

In each of the cases we have checked, the PH-DMRG leads to a significant 
improvement of the ground-state energy approximation from that of the HF 
method. For example, while the HF results, for a 4 x 6 lattice occupied by 10 
electrons with NN interactions, show a discrepancy of 4.5 percents, the PH- 
DMRG results, which were obtained by keeping up to p = 100 block-states 
between successive iterations, decrease it to ~ 1.5 percent from the exact 
solution. We have compared this method to two other methods which were 
reported in Ref. [ 116 ], and which can be used to approximate the ground- 
state energy with a truncated Hilbert space. The obtained PH-DMRG error 
rate is better than the method we've denoted as energy-cut, and is similar 
to that of the generation- cut technique. 
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Nevertheless, for a full comparison one must also consider the feasibility 
of these methods when a treatment of larger systems is required. To under- 
stand the difference between the methods, let's start with the small lattice 
size we've treated, 4x6. The energy-cut, based on 18 energy levels, and 
the generation- cut, with up to 3 particle- hole generations, require the diag- 
onalization of matrices of sizes 43, 758 and 47, 916, respectively. The largest 
superblock diagonalization in the PH-DMRG process, for p = 120, was for a 
matrix size of 13, 494. 

When a larger lattice is treated, and in order to get the same accuracy 
as for the small system, the energy-cut and generation-cut methods must 
include more levels and generations. Moreover, even was it sufficient to take 
the same number of levels and generations as in the small case, the size of 
matrices would have grown exponentially with the lattice size. Therefore, the 
lattice enlargement makes these methods infeasible even for modest lattice 
sizes. In the PH-DMRG method, on the other hand, the size of the matrix 
may remain constant, since it depends on the block size, and not on the 
system size. 

Yet, for larger systems, the number of single-particle states is larger, 
and the discrepancy of the PH-DMRG is expected to increase, unless more 
block states are constantly kept. Therefore p, and thus the matrix size being 
diagonalized, are expected to increase in the PH-DMRG method as well. 

The largest matrix size needed to be diagonalized in the PH-DMRG pro- 
cess is shown in Fig. 15.71 as a function of p. The dependence of Mmax on p, 
for large values of p, was empirically found to be linear [61] . From our results 
one can see that although a linear fit of the large-j9 points (orange line) is 
possible, yet a power law including the entire p regime (black line) seems 
more appropriate, resulting in M^ax ~ p^-^^. In any case, it is clear that the 
matrix size is less than quadratic in p. Furthermore, since the largest matrix 
size used in our current PH-DMRG application is still much below the tech- 
nology limit, its increase should not be a problem. It is thus clear that in 
principle the PH-DMRG is capable of treating larger systems. Indeed, initial 
studies we have already performed show that the PH-DMRG is feasible for 
systems of the order of 10 x 10. 



* * * 
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Figure 5.7: The maximal size of the superblock Hamiltonian needed to 
be diagonaUzed, for the 4x6 lattice with 10 electrons. The orange line is 
a linear fit of the p > 80 points, while the black line is a power law fit in 
the entire range. 



Beyond the ground-state energies, we have used the improved PH-DMRG 
method in order to calculate the addition spectrum, and we have checked its 
accuracy. The best accuracy we have obtained is for (A2), in which a rapid 
convergence is exhibited, and a very small number of block-states is sufficient, 
for both short-range and long-range interactions. The accuracy of the 6A2 
calculation was pretty poor in the NN case, resulting in more than 5% error 
and very slow convergence for p = 60. Even though, the PH-DMRG results 
of 6A2 are much more accurate than the corresponding HF results. Moreover, 
for the case of long-range interactions, the 6A2 calculation was much better, 
getting an accuracy of ~ 0.8 percent. The A2 discrepancy, averaged over the 
different realizations, has obtained in the long-range case ~ 0.13 percents of 
error, much better than the ~ 1.7 percents of the HF results. 

A comprehensive study of the addition spectrum was performed a few 
years ago, using the self-consistent HF method, on 2D lattices of the order 
of 10 X 10 |118] . It was found that for strong interactions the fluctuations 
of A2 do not scale with the mean level spacing A for both short-range and 
long-range interactions, as opposed to a single-parameter scaling argument. 
In the case of short-range interactions, a saturation of 6A2 was obtained 
for very high values of interaction and it was explained as an appearance of 
charge density modulations. For Coulombic interactions the increase of 6A2 
was shown to be faster than linear in A. 

Based on our results, however, one might question the ability of the HF 
method to provide an accurate approximation of the addition spectrum. Our 
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results point out that the usage of the HF method may lead to quite high 
error rates, especially for short-range interactions. It can thus be interesting 
to investigate the addition spectrum for such systems using the PH-DMRG, 
which is appropriate for such system sizes, and which can improve the ac- 
curacy significantly A PH-DMRG research of the addition spectrum is thus 
left for a future study. 



Chapter 6 

Two-Electron Magnetization in 
Quantum Dots 

In this chapter we study the magnetization of an interacting quantum dot 
occupied by spin 1/2 electrons, in the presence of spin-orbit couphng. Mod- 
ehng the dot by a two-dimensional lattice, we utilize exact calculations for 
lattices with a small number of electrons. For an A^-sites lattice, which is oc- 
cupied by rif, electrons, the Hilbert space size is M = pj^), so that if the size 
of the system is not too large, and it is occupied by relatively small number 
of electrons, one can still use an exact diagonalization method and get the 
exact lowest wave-functions, even when interactions between the electrons 
are considered. 

We calculate the expectation value of the spin operators in the ground- 
state, and by a comparison of the system's energy with and without a mag- 
netic field, we also calculate the g-factor. For the case in which there are 
two electrons occupying an interacting quantum dot, we find a level crossing 
between the two lowest many-body eigenfunctions as a function of the spin- 
orbit scattering rate, resulting in a finite magnetization of the ground-state. 
This is a clear evidence for the importance of the interplay between spin-orbit 
scattering and interactions, and can have a significant influence on g-factor 
measurements. 

6.1 Introduction 

The effects of spin-orbit (SO) coupling on the energy spectrum of quan- 
tum dots (QDs) have attracted notable attention in the recent years [1191 
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\120\ 11211 1122] . Much work has been concentrated, both experimentally and 
theoretically, on the magnetization of mesoscopic samples. Specifically many 
studies focus on g-factor measurements, and try to understand them theoret- 
ically. The g-factor is defined through the splitting of the Kramers' doublets 
|123l [To] in the presence of a weak magnetic field, e\^^ — ef'^ = ±^gfiBH, 
where e[^^ i^f^) is the i-th single-particle energy level with (without) a mag- 
netic field, /iB is the Bohr magneton and H is the magnetic field. 

For free electrons the g-factor is constant, g = 2, and this value is more or 
less correct also for bulk measurements in various metals ^119j . However, in 
experiments done on mesoscopic samples, values which are significantly less 
than the free value of the g-factor were obtained. Moreover, large fiuctua- 
tions in the measured values were seen. These findings attracted theoretical 
attention, and resulted in two analytical studies which have obtained, within 
the framework of the random matrix theory (RMT), a description of the g- 
factor probability distribution in the presence of SO coupling and disorder 
but without interactions |124^ I125j . In a more recent study, the statistical 
properties of these distribution functions were related to several physical ob- 
servables |126j . According to these studies, the SO coupling infiuences the 
probability distribution of the g-factors of the discrete energy levels. The dis- 
tribution function was shown to be universal, where the width is expressed 
in terms of various physical parameters. The presence of strong SO coupling 
and disorder causes the g-factor to fiuctuate from sample to sample, even 
when both have the same characteristic strength of SO coupling. Moreover, 
the g-factor is expected to fiuctuate also between different levels of a specific 
sample, according to the RMT distribution function. 

The definition of the g-factor through the splitting of the Kramers' dou- 
blets thus relates the g-factor to the change in the energy of a specific energy 
level when a weak magnetic field is applied. Therefore, 
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where Qi^a is the g-factor of the i-th level, with spin cr. The spin index 
(7 G {+, — } is used to relate the states \i, +) and \i, — ) via a time reversal 
operation. In the absence of a magnetic field each level is two-fold degenerate 
(Kramers' degeneracy |123l fT0]). and this degeneracy is lifted by the magnetic 
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field, which causes the energy of one of the levels to increase and the energy 
of the other to decrease. Therefore, Qi^^ as defined by this formula can have 
either sign, depending on the direction of the energy change. The ground- 
state energy always decreases when a magnetic field is applied, thus the 
g-factor of the ground-state obtained by Eq. (16.11) is positive. Usually, the 
value of g does not depend on the spin index, at least to zeroth order in H, so 
that one can denote the g-factor of the i-th level as ±gi, with the convention 
that gi > 0. 

According to Eq. (16.11) , a measurement of the g-factor should compare the 
specific energy level before and after the magnetic field is applied. However, 
practical measurements are usually related to the total energy of the system, 
and not to that of a specific level. Nevertheless, if there is no interaction be- 
tween particles, the change of the total ground-state energy due to magnetic 
field may be equivalent to that of the highest filled level. This is the case 
when the number of electrons rig is odd, i.e. Ue = 2p + 1, with an integer 
p. The lowest rig — 1 levels are composed of p Kramers' pairs, where in each 
pair one level increases and the other decreases in the presence of a magnetic 
field, so that their total contribution vanishes. Therefore, the motion of the 
highest level, p+1, which is singly occupied, will determine the g-factor, so 
that one can write 

gine) = 77 ^, (6.2) 

where g{ne) denotes the g-factor of the ground-state with rig electrons, and 
Egf\ne) represents its energy in the presence of a magnetic field H. 

In addition, when the number of electrons is even, the total ground- 
state energy is not expected to change when a magnetic field is applied, since 
all the filled levels divide into pairs, in which the movement of one level is 
compensated by the othei0. Therefore, for an even number of electrons, a 
calculation of the g-factor using Eq. (16.21) gives g = 0. One can thus use 
Eq. (Q practical definition of the g-factor. 



^Actually, in the presence of SO coupling and magnetic field, there are some second 
order corrections to the energy, which lead to first order corrections to the g-factor. This 
issue will be addressed in section l673l However, to first order in H the ground-state energy 
of an even-occupied non-interacting system is constant, and g — 0. 
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However, returning to the experimental point of view, one should notice 
that when the measurement of the energy is indirect, an interpolation of the 
results is also required. For example, using tunneling spectroscopy one can 
measure the position in which a conductance peak of a QD occurs. Such an 
event relates the energies of the QD with Ug — l and electrons, with the gate 
voltage Vg, by the relation (see section [1.1.21) eVg = Egg^rig) — Egsirif. — 1). 
When a magnetic field is applied, the peak will move as a function of H. 
Therefore, by denoting the measured g-factor by one can analyze the peak 
motion in order to determine the g-factor, by calculating 

2 [eVg{Q) - eVgiH)] 
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(6.3) 
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Since either He or rig — 1 is even, its corresponding g-factor vanishes, and 
thus g is equivalent to the other. Namely, g = g{ne) or g = —g{n(. — 1). 
Actually, since each peak is divided, in the presence of a magnetic field, into 
two peaks moving in opposite directions, an extraction of g from successive 
peaks results in a set of the single-particle g-factors, i.e. gi, —gi,g2, —92, ■ ■ ■ ■ 
Indeed, tunneling spectroscopy measurements have obtained many rea- 
sonable results for g-factors of nano-particles. For example, several experi- 
mental studies of metallic three dimensional nano-particles have shown the 
reduction of the measured g-factor as a function of the spin-orbit scattering 
rate, in accordance with the RMT predictions. For Aluminum nano-particles, 
in which the SO coupling is negligible, the measured g-factor values are 
approximately those of free electrons {g = 2) |121j . while for Gold nano- 
particles, in which the SO coupling is strong, the measured g-factors were in 
the range of 0.28 — 0.45 [122]. Furthermore, by extracting several g-factors 
from each sample, Petta and Ralph have succeeded to present an impressive 
confirmation of the theoretical RMT distribution function [127] . Neverthe- 
less, it should be mentioned that while the average g-factor is expected to be 
reduced by the SO coupling for three dimensional samples, the RMT predicts 
an enhancement of the g-factor in two dimensional (2D) samples [125] as a 
function of SO coupling. Yet, no experiment which measured this increase 
has been performed to date. 
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However, as we have mentioned, Eq. (16.21) is entirely equivalent to Eq. (16. ip 
only when the system is non-interacting. Once interactions between the elec- 
trons are important, it should be emphasized that Eq. (16.21) is a definition 
of a many-particle g-factor, which depends on the total magnetization 
of the ground-state wave-function. For example, one can obtain values of 
g-factors which are larger than 2, even in the absence of orbital effects (e.g., 
when the magnetic field is in-plane), a phenomena that cannot happen for a 
single-particle g-factor. 

Nevertheless, the RMT-based theoretical studies cited above were per- 
formed in the context of one-body levels, neglecting any effect of the electron- 
electron interactions. Indeed, by adding to the RMT Hamiltonian an inter- 
action term, using the constant interaction model, an increase of the g-factor 
fluctuations was reported |128[ 1129] . It was shown that the interactions result 
in a possibility of getting non-trivial spin values in the ground-state, and ac- 
cordingly in an optional enhancement of the g-factor, to values greater than 
2. 

Although the theoretical studies of Refs. |128[ 1129] were done for an odd- 
electron occupation, their results suggest the possibility of a non-trivial spin 
polarization for the even-electron well. Similar phenomenon was 

found in disordered dots without spin-orbit coupling and with inflnitely large 
Hubbard interaction, where occupation of an even number of electrons caused 
non- vanishing spin values |13U] . If, for any reason, this is the case, and the 
g-factor of an even-electron ground-state differs from zero, then the quantity 
measured in tunneling spectroscopy may not equal the single-level g-factor 
nor the many-particle g-factor. In such it should be related to the 

difference between two many-particle g-factors, as shown in Eq. (16.31) 

In this chapter we investigate the ground-state magnetization properties, 
such as the spin polarization and the g-factor of QDs with an even number of 
particles. Without interactions, such states have g = 0, a.s well as {Sz) = 0, 
as predicted. Nevertheless, we show that the interplay between spin-orbit 
scattering and electron-electron interactions may result in a level crossing 
(LC) between the two lowest many-body levels. When these states are close 
in energy, the magnetic fleld splits them into two polarized states with a flnite 
magnetization. Therefore, there exist a possibility to have non-vanishing (S^) 
and g-factor in the two-particle ground-state. 
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6.1.1 Chapter's Outline 



The rest of the chapter is organized as follows. In the next section we de- 
scribe the model Hamiltonian we use in order to incorporate, beside the 
magnetic field, both SO coupling and interactions between electrons. In sec- 
tion 16.31 we present results for a non-interacting system, for both single- and 
double-occupation, which are shown to reproduce some known ground-state 
properties. We find that there are specific values of the SO coupling strength, 
in which the Kramers' doublet remains degenerate even when a magnetic field 
is applied. In such points both states of the doublet have (Sz) = 0. 

The results for the case in which interactions are also considered are 
presented in section 16.41 Our results point out that a finite magnetization 
can be obtained for systems with an even-particle occupancy. In section 
16.51 we discuss the experimental relevance of this finding, and we show that 
it might affect g-factor measurements. In the last section we conclude and 
address some future possibilities to continue the research. 



In order to model the QD we choose a tight-binding description of a finite 
2D lattice with A columns and B rows (the number of sites is denoted by 
N = AB), with open boundary conditions, which is occupied by Uf, spin 1/2 
electrons. In some sections of the chapter, only point interactions are con- 
sidered, i.e., the Hubbard term which couples up and down spins occupying 
the same lattice site. In other parts we add to the Hubbard term either 
nearest neighbor (NN) or Coulomb interactions. In addition, spin-flips dur- 
ing hopping processes are possible, with a finite probability, as a result of a 
coupling between the spin degree of freedom and the orbital motion. Sep- 
arating the interactions from the free part, one can write the Hamiltonian 
as i^QD = Hq + Hint, where the free part can be divided, if disorder effects 
are neglected, to a hopping term and a Zeeman term, i.e. Hq = Hhop + Hb- 
Each of these terms was discussed in details in section II. 1.11 The hopping 
part of the Hamiltonian is thus 
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where aj^na {(^m,n,a) is a creation (annihilation) operator of an electron with 
spin a in the lattice site placed in row m and column n. The matrices 
and Vy are defined (in the absence of a magnetic field) by |12] 

where Vi (V2) is the hopping matrix element, for events which conserve (fiip) 
the spin. 

When a perpendicular magnetic field is applied, it adds a phase to the 
hopping matrix element. We take the field direction to be perpendicular 
to our 2D sample, i.e. along the z axis, and we choose a gauge in which 
the vector potential is A = Hyx. With that gauge, one has to modify the 
hopping elements is the x direction (see section [1.1.11) . so that — > Kce*^*", 
with m being the row number and 9 = ^^^^^ . Here s is the lattice constant 
and 00 = hc/e is the magnetic fiux quantum. In both cases, i.e. with and 
without a magnetic field, the overall hopping amplitude, t = \/Vi + V2, is 
taken as the energy scale of the problem. In other words, all energy terms 
are expressed in units of t. 

The strength of the SO coupling can be expressed by the ratio between 
the the spin-fiip amplitude and the total hopping element (excluding phases). 
Using a dimensionless parameter A = , = Vo/t, we examine the entire 
range of A, between very weak (A 0) and very strong (A < 1) spin-orbit 
coupling. A similar approach is usually utilized within the RMT framework, 
by writing H = {1 — a)HGOE + oiHgse- Here GOE (GSE) denotes the 
Gaussian orthogonal (symplectic) ensemble, which corresponds to the case 
in which the SO coupling is very weak (strong). Changing a from to 1 
modifies the Hamiltonian between these two limitj^. 

With the choice of a perpendicular magnetic field, the Zeeman term in 
the Hamiltonian is diagonal, and can be written as 

Hb = I^bH ^ f^a]n,n,a«m,n,a, (6.6) 
m,n,cr 

where a = ±1. 

^Note that the definitions of A and a are different. A strong SO coupUng which is related 
to the symplectic ensemble of RMT (a — > 1), corresponds, in our model, to A w 0.5 — 0.7. 
The physical relevance of higher values of A is questionable. 
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The Zeeman energy can be related to the hopping phase 6 and to the 
hopping amphtude t by the following consideration. The phase ^ is a dimen- 
sionless parameter, measuring the magnetic flux throughout a lattice unit 
cell, in units of the quantum flux 0o- One can express the absolute value of 
the Zeeman energy as fisH = /iB0o2^- Substituting the physical constants 
/^s0o = n^o being the electron mass, and using the relation t = ^m^is^ ' 
where meg is the effective mass, one gets /XB-ff = = For the 

metallic nano-particles used in several experiments |121l 11221 1127j . the ratio 
between the effective mass and the electron mass is close to unity. We will 
assume meg = mo; Deviations from this value will not affect our main re- 
sults. Finally, since all energies are measured in units of t, the strength of 
the magnetic field ^sH/t determines exactly the hopping phase. 

The Hubbard term results in 

m,n 

where Uh is the Hubbard interaction strength. When an increase of the 
range of interactions is considered, by either NN or Coulomb interactions, 
one of the following terms is added to the Hubbard term: 

<mi,ni;m2,n2>;(Ti,cr2 

C = E — - — 

mi ,ni 7^m2 ,n2 ;fTi ,(T2 



'"mi,rai '^m2,n2 



'^mi ,ni ,(Ti '^m2 ,n2 ,(72 '^"12 ,«2 ,cr2 '^'mi ,ni ,cri ; 



where < mi, rii; m2, '^2 > denotes that the sites mi,ni and m2,n2 are NNs, 
and [/nn is the interaction strength between NN. For the Coulomb inter- 
actions, the distance between sites mi,ni and m2,n2, expressed in lattice 
constant units, is denoted by \rmi,ni — rm2,n2\, and Uc is the Coulomb inter- 
action strength between sites which are one lattice constant apart. 

In order to calculate the spin polarization of the QD we apply a weak 
magnetic field along the z axis and calculate (5*^) for the lowest levels. For g- 
factor calculations, we compare the ground-state energies with and without 
the magnetic field for each sample, and use Eq. (16.21) . We use fisH/t ~ 
10~^ — 10"'^, for dots in which the mean level spacing is of the order of O.lt. 
For an experimental system in which the mean level spacing is 0.1 — 1 meV, 
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it is equivalent to a magnetic field of 10 — 1000 G, in correspondence with 
practical measurements |127j . 

The Hamiltonian Hqb is exactly diagonalized using Lanczos procedure, 
for lattices of up to 11 x 10 sites, occupied by 1 or 2 electrons. A discussion 
of the lattice sizes will be given in section 16.51 It will be shown that some 
of the significant results, to be detailed in the following sections, are not an 
artifact of small sizes. From checking the size dependence of these results 
one can see that they will not disappear, even for larger QD sizes. Moreover, 
they are expected to be even more pronounced. 

6.3 Non-Interacting Electrons 

We start by presenting the non-interacting results for the spin polarization, 
by taking Uh = f^NN = Uc = . The expectation value of the operator 5*2 
in the low lying states, i.e. (Sz), is calculated when a weak magnetic field 
{fiH <^ t) is applied along the z axis. 

The results for {S^) as a function of the SO parameter A, of singly- 
occupied states in a 8 X 7 lattice, are shown in Fig. 16. 1[ When it is needed, we 
denote by {Si"^^) the expectation value of the operator Sz in the m-th eigen- 
function (m = 1 being the ground-state). In places where expressions with 
different numbers of particles are related, we use {Si^\ne = n)) to denote 
that the expectation value is evaluated in the n-particle Hilbert space. 

We start with the single-particle levels, taking rie = 1. Without the mag- 
netic field, all single-particle states (and in particular the ground-state) are 
doubly-degenerate (the Kramers' degeneracy) [123^ [TU] . When a magnetic 
field is applied, it splits this degeneracy, and one gets to zeroth order in the 
magnetic field, {S^z"*) = —{Si^^). For A ^ 0, \{Sz)\ |, and increasing the 
SO coupling leads to a decrease of K^^)!, as can be expected. For different 
levels one gets similar, although not identical curves, with the same qualita- 
tive limits for weak and strong spin-orbit coupling. Results for the g-factor 
calculation are similar, i.e. for A ^ we get (for the ground-state) g = 2, 
and as the SO strength increases g drops monotonically towards g = 0. 

In the insets of Fig. 16.11 we zoom into the regime of strong SO coupling, 
showing the difference in the energies of the two lowest states (upper inset), 
and each of their (Sz) (lower inset). Both plots point towards LCs between 
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Figure 6.1: The spin projection (Sz) of the lowest two single-particle 
levels are shown as a function of the spin-orbit coupling strength A, for a 
non-interacting system of 8 x 7 sites. Insets: zoom into the strong spin-orbit 
regime shows signs of level crossings. The upper inset shows the energy 
difference between the first two levels when a magnetic field is applied 
(notice the logarithmic scale), and the lower one shows {Sz) of each. 



the first two single-particle levels. At the crossing points both states have 
{Sz) = 0, so that the magnetic field does not break Kramers' degeneracy, 
an unusual result. It should be noted, however, that the crossing points are 
obtained for values of A > which correspond to the situation in which a 
spin-flip hopping process is more probable than a spin-conserving one, and 
the existence of such a regime in practice is improbable. 

For n = 2, without interactions, the many-body state is a Slater de- 
terminant of the single-particle states, and since the operator Sz is addi- 
tive, one can write {Si^\ne = 2)) = {Si^\n(, = 1)) + {Si^\ne = 1)), and 
{Si^\ne = 2)) = {Si^\ne = 1)) + {Si^\ne = 1)). To zeroth order in the mag- 
netic field, as noted above, the two contributions to {Si^\ne = 2)) cancel 
each other, and this term vanishes. The results are shown in Fig. 16. 2[ The 
results for the g-factor calculation are similar, i.e. for the ground-state we 
get g ^ 0, while for the first excited state (7 = 4 for A = 0, and it decreases 
monotonically towards g = when A increases. 

Since for the g-factor measurements the limit H ^ is taken, the ap- 
proximation {Si^\ne = 2)) ^ is usually sufficient. However, as can be seen 
in the inset of Fig. 16. 2[ there is also a higher order term. For weak spin-orbit 
coupling and weak magnetic field, using first-order perturbation theory in 
both A and H, it was shown that there is an additional contribution to Sz, 
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which is proportional to X^H |131] . This explains the quadratic increase of 



{S. 



(1), 



rir. 



2))/^bH as a function of A which is shown in Fig. I6.2( inset) for 



weak SO coupling, in which the perturbation theory is valid. For moderate 
values of SO coupling, our results suggest that this trend is reversed, and 



(1), 



Tip 



2)) / hbH starts to decrease. 




Figure 6.2: The spin projection {Sz) of the lowest two doubly-occupied 
levels are shown as a function of the spin-orbit coupling strength A, for a 
non-interacting system of 8 x 7 sites. Inset: the spin polarization of the two- 
particle ground-state, which has a first-order dependence on the magnetic 
field H. For weak SO coupling, it also has a quadratic dependence on A 
(the dashed brown line represents a quadratic fit for small A). 



Returning to the points of LC shown in the insets of Fig. 16.11 it is impor- 
tant to notice that such crossings occur between states which are the time 
reversal of each other, or, in other words, states which belong to the same 
Kramers' pair. No such crossings occur between states which originate from 
different pairs, although the single-particle level-spacing is significantly re- 
duced by the SO coupling. Yet, the energy difference between subsequent 
Kramers' pairs is much larger than the energy contribution of the weak mag- 
netic field we apply. Correspondingly, for the double-occupation case, the 
lowest two states (those which are presented in Fig. 16.20 should not exhibit 
any crossing. This can change once interactions are considered, as will be 
shown in the next section. 

Whereas the exact point of the maximum in {sP(2))/fiBH is different for different 
iattice sizes, the quaiitative shape for ali the sizes checked was the same. 
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6.4 Interplay Between Spin-Orbit Coupling 
and Interactions 

We now move to study the effect of interactions, by turning on the Hubbard 
interaction term, for which an interaction energy Uh > is paid for a couple 
of electrons occupying the same lattice site. Calculating the ground-state 
energies of the two lowest doubly-occupied states, one finds that there is a LC 
between these states, a feature which does not exist for the non-interacting 
case. In the non-interacting case different levels may approach each other 
when the SO coupling increases, yet the minimal distance between them is 
much larger than the magnetic energy. The presence of interactions enhances 
this tendency, towards the situation in which a LC is possible. This crossing 
happens at a certain value of the SO coupling, i.e. at A = Ac, and in its 
vicinity, the expectation value of S"^ switches between these states, as can be 
seen in Fig. 16. 3[ 
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Figure 6.3: Typical results of the level crossing of the two lowest doubly- 
occupied states. The results shown were obtained for a system of 8 x 7 
sites, with Uu = 10. Upper panel: the energy difference E2 — Ei (without 
a magnetic field) is shown as a function of the spin-orbit coupling strength 
A (notice the semi-log scale). The dip shows the crossing point. Lower 
panel: the switch of (5"^) (in the presence of a magnetic field) between 
these two states, which occurs at the same place. Note the tiny scale of A. 

Looking at the energy curves and the switching of {S"^)-, one would naively 
expect that the magnetization properties, e.g. {Sz) and the g-factor, will 
switch as well at Ac. However, as the energies of these two states become 
close enough to each other, the energy associated with the magnetic field 
becomes more and more important, resulting in a polarization of the spins 
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of both states. This leads to an enhancement of (Sz) and the g-factor values 
in the crossing region. As can be seen in Fig. 16.41 both (5*2) and the g-factor 
can reach significant values. 

0.06 P ' \ ' \ ' I ' ~ 




I > \ I \ > I > I 

0.676 0.677 0.678 0.679 0.68 

X 



Figure 6.4: Typical results of the spin polarization (Sz) (upper panel) 
and the g-factor, calculated using Eq. (j6.2p (lower panel), of the two lowest 
doubly-occupied states, in the regime of the level crossing between them. 
The results shown were obtained for a system of 8 x 7 sites, with Uh = Wt 
and ^ibH = 10~^t. Note the tiny scale of A. 

Looking at the particle distribution throughout the lattice can shed some 
light over the mechanism of the LC. For most values of A the first two doubly- 
occupied states are different in their density distribution, as can be seen in 
Fig. 16.51 However, in the vicinity of the LC, these states are identical in their 
spatial components (Fig. 16. 6p . Yet, their spin degree of freedom gives rise to 
a polarization of both states, in opposite spin directions (Fig. 16. 7[) . One may 
conclude that a specific combination of the SO and the interaction strengths 
can lead to a separation of the spatial and the spin degrees of freedom. In 
such the lowest two states are identical in their spatial coordinates, 

whereas the spin degree of freedom is responsible for their polarization. In 
other words, states which experience a LC are time-reversal of each other, 
similarly to the non-interacting case. 

We thus see that in the vicinity of a LC a finite value of the g-factor can 
be obtained, in contrast to the ordinary assumption of g = for a doubly- 
occupied system. In general, the vanishing of the g-factor results from the 
quadratic dependence of the ground-state energy on the magnetic field. Such 
a dependence is shown in the upper panel of Fig. 16.81 for an arbitrary value 
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Figure 6.5: The electron distribution of the first two doubly occupied 
states (left - ground state, right - first excited state) for a 8 x 7 lattice with 
Hubbard interactions of strength Uh = lOt and SO coupling A = 0.550. 




Figure 6.6: The same as Fig. 16. 5( with A = Ac = 0.678, resulting in 
identical charge distributions. 




Figure 6.7: The electron distribution of the ground state of Fig. 16.61 

separated between spin-up (left) and spin-down (right) electrons. In the 
first excited state these distributions are exchanged. 
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of the SO coupling. However, near a LC point, as shown above, the ground- 
state has a finite spin polarization. As a result, the dependence of the energy 
on the magnetic field is linear (lower panel of Fig. 16.81) . and the g-factor is 
finite (lower panel of Fig. 16. 4p . The clear linear dependence in the exact 
point Ac is actually limited to the region in which (IbH is greater than the 
energy difference of the two many-particle states. The same restriction holds 
for the peaks in (Sz) and g, which are thus getting wider as the magnetic 
field is enhanced. 




Figure 6.8: The dependence of the energy on the magnetic field is com- 
pared between the regime of a level crossing (lower panel), to another ar- 
bitrary point (upper panel). The results shown by symbols were obtained 
for a system of 8 x 7 sites, with Uh = 10, and the solid lines represent 
quadratic (upper panel) and linear (lower panel) fits. 

From these results one can conclude that whereas the g-factor vanishes by 
definition for most values of A, it has a finite value near Ac- Calculating the 
g-factor using Eq. ( 16. 2p leads to the typical results shown in the lower panel 
of Fig. 16. 4[ As will be discussed in the next section, such g-factor values can 
be significantly large, and thus they cannot be neglected. 

6.5 Experimental Relevance 

The enhancement of the g-factor discussed in the previous section, was ob- 
tained for various system sizes. However, since the exact diagonalization 
technique used is limited by size, the scalability question, i.e. the question 
whether such a finite g-factor can be experimentally measured, is important. 
The peak, in both (S^) and the g-factor, can be characterized mainly by 
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two properties, namely the peak height and its width. In addition, the peak 
occurrence is characterized by the value of Ac, and by the corresponding 
interaction strength and lattice size for which it occurs. 

To answer the scalability question one should check if for realistic sample 
sizes with a reasonable strength of interaction the crossing point Ac is small 
enough to be obtained by realistic doping with magnetic impurities. In ad- 
dition, one has to find whether the width and the height of the predicted 
g-factor peak at this point are experimentally measurable. 
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Figure 6.9: The g-factor of the doubly-occupied ground-state as a func- 
tion of the SO coupling strength A, for lattices sizes of 8 x 7 (black), 9x8 
(red), 10 X 9 (green) and 11 x 10 (blue). Different panels correspond to 
different type of interactions: Hubbard interactions with Uh = Wt (upper 
panel), Hubbard and NN interactions with Uh = lOi and Unn = 5t (mid- 
dle panel), Hubbard and Coulomb interactions with Uh = lOt and Uc = 5t 
(lower panel). 



Since a substantial enlargement of the system is not numerically possible, 
in the following we present results for different system sizes, from which 
the trend can be clearly understood. For example, in the upper panel of 
Fig. 16.91 the g-factor peak is shown for system sizes ranging from 8 x 7 to 
11 X 10. As can be seen, the value of Ac decreases with increasing system 
size, suggesting that for a sufficiently large system size the crossing occurs 
for moderate value of the SO coupling. On the other hand, since the peak 



6.5. EXPERIMENTAL RELEVANCE 



137 



height decreases slowly with increasing system size, the question if it will not 
be negligible for realistic system sizes remains open. 

Nevertheless, the effect of other types of interactions can change this 
picture. Using NN interactions (in addition to the Hubbard term) results in 
much higher peaks, so that there is a larger chance to find a finite g-factor 
for larger systems. This is shown in the middle panel of Fig. 16.91 In the lower 
panel we show the results obtained when Coulomb interactions were added 
to the Hubbard term. Although the peak heights are smaller, they become 
much wider, and most important, the peak height increases with system size. 

Further investigation of the Coulomb case shows that the maxima in the 
g-factor curves are not accompanied by LCs between the first two levels. Yet, 
they occur when the energies of these levels are close enough. Furthermore, 
these maxima, and their different shape, result from the second order term 
(i.e., the term) in the energy. The corresponding g-factors are thus of 
first order in the magnetic field (and not constant as in the Hubbard and 
NN cases). Nevertheless, significant values are obtained, even for a magnetic 
field as weak as fiH/t ^ 10~'^, corresponding to 10 — 100 Gauss. 

Moreover, since the maxima are the result of avoided crossings, and the 
energy difference between the levels is expected to reduce when larger samples 
are treated, the values of the g-factor in the maxima region should further 
increase. This is in contrast to the other cases (Hubbard, NN) where the peak 
results from a true LC and decreases with increasing system size, probably 
due to a decrease in the inter-level Hamiltonian matrix element caused by 
the applied magnetic field. 

In order to get insight into the other differences between the three cases, 
namely the LC and the sharp g-factor peak in the Hubbard and in the NN 
cases, which are absent in the Coulombic case, the role of the spin component 
is further explored. By considering f/|| = Un = Ui and U-^i = f/j^f = U2 (for 
either Unn or Uc), one can check whether these phenomena are observed 
for different parameter regimes, and specifically for the limits of parallel spin 
interactions {U2 — > with a finite Ui), and anti-parallel spin interactions 
(f/i and U2 is finite). 

For the NN case, we find that the g-factor peak is enhanced when the 
interactions are only between anti-parallel spins, while it disappears for the 
case in which only interactions between parallel spins are considered. For 
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the Coulombic case, when only anti-parallel interactions are used, a sharp g- 
factor peak does appear. However, when parallel interactions are considered, 
either with or without the anti-parallel ones, the sharp peak disappears. In 
all cases an appearance of the sharp g-factor peak is accompanied by a LC. 
The different cases are summarized in Table 16.11 





Spin independent 


Parallel 


Anti-parallel 




(f/i = U2) 


{U2 ^ 0) 


(f/i - 0) 


NN 


LC 




LC 


C 






LC 



Table 6.1: Level crossing occurrence for different interaction types 



These results can be understood in the following way. As we have previ- 
ously shown, the crossing levels are the time reversal of each other, and the 
spin degree of freedom is responsible for their spin polarization. Therefore 
the question whether a LC can occur is crucially related to the possibility 
to polarize the lowest states. The polarization of these states might become 
improbable when interactions between parallel spins are considered. 

Therefore, when there is only anti-parallel interactions, for both inter- 
action types (NN, Coulomb) a polarization of the ground state is possible, 
and thus a LC does occur. In the opposite case, when the interactions are 
only between parallel spins, the probability of spin polarization is reduced 
because of the interaction. For the case of Coulomb interactions, the polar- 
ization of the lowest states is totally blocked, and a LC cannot occur. In the 
NN case, however, since the interaction is only short ranged, single-particle 
states which are spatially separated can be combined, in a rough approx- 
imation, to a polarized many-particle state. However, a two-particle state 
composed of anti-parallel spins is energetically preferable, for any strength 
of A, so that a LC of the lowest states does not occur. 

In the case of spin-independent interactions the difference between Coulomb 
and NN interactions stems from the interaction range. The presence of 
both parallel and anti-parallel interactions types causes the re-appearance 
of a polarized ground-state in the NN case, combined of spatially separated 
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single-particle states. On the other hand, the long range of the Coulomb 
interactions prevents the ground-state polarization, and thus a LC between 
the lowest two levels does not occur. It should be noted, however, that a 
strong avoided crossing does appear (lower panel of Fig. 16.91) . In addition, 
a polarization of an excited state, as well as LC between excited states, are 
still possible. 

6.6 Conclusions and Future Prospects 

In this chapter we have shown that the combination of interactions and spin- 
orbit scattering can cause unexpected magnetization of states with an even 
number of electrons. We have also shown that for realistic sizes of QDs, 
such a result can be experimentally observed, and might be relevant for 
understanding some measurements. 

As we have noted, one of the popular methods for measuring the g-factor 
is by tunneling spectroscopy. The number of electrons in the dot in such 
an experiment is changed by one during each tunneling event, involving a 
transition between an even electron number and an odd one. According to 
our results, such a measurement might present the result for the difference 
in g-factor between the two states, Eq. fl6.3p . If the even-electron state has a 
non-vanishing g-factor, like in the vicinity of the LCs we have presented, the 
measured quantity g may not equal the g-factor of the odd electron state, to 
which it is usually attributed. 

In such trace of the LC may be seen experimentally. In the 

regular case (as opposed to the LC scenario), the two levels which belong to 
the same Kramers' doublet have the same g-factor, and the motion of the 
two energies as a function of a magnetic field is symmetric. However, in the 
region of a LC, the two levels get contributions from different even-particle 
states. Explicitly, the p Kramers' pair is divided by the magnetic field to the 
measured values g{2p — 1) — g{2p — 2) and g{2p) — g{2p — 1). In general, this 
motion, as a function of the magnetic field, is not symmetric. An example is 
presented in Fig. I6.10[ As one can see, the most clear non-symmetric motion 
is obtained for A ~ Ac (left panel), but such a motion can be seen for a 
region in its vicinity as well (middle panel). Far enough from this region 
(right panel) the symmetric motion reappears. 
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Figure 6.10: The motion of the first Coulomb peaks (rig = 1,2) as 
a function of the magnetic field is shown for a lattice of 11 x 10 with 
Uh = Wt and Uc = 5t (for which ^ 0.31). 

* * * 

As a last remark we note that the model we have used in this chapter has 
neglected any kind of disorder inside the QD. This is in contrast with the 
usual experimental configuration, in which the fabrication of clean samples is 
difficult. It is thus interesting to check the influence of disorder on the two- 
particle g-factor peak described above, for example to investigate how its 
place and shape vary between different realizations. We leave this question 
to future research. 



Chapter 7 
Summary 



In this thesis we have examined two types of mesoscopic systems which are 
commonly used nowadays in innovative physical studies. The first type 
we have investigated is a quantum dot (QD) which is coupled to a one- 
dimensional (ID) wire, and the second type is an isolated two-dimensional 
(2D) QD. 

When a QD is coupled to one end of a semi-infinite ID lead, it may have an 
important influence on the wire's characteristics. Using the numerical density 
matrix renormalization group (DMRG) method we have shown in chapter 
[3] that the thermodynamic properties of the entire system are sensitive to 
the dot properties: its energy level and the strength of its coupling to the 
wire. The wire itself can be in one of three phases, the ferromagnetic (FM), 
the Tomonaga-Luttinger liquid (TLL) and the charge density wave (CDW) 
phases. When the wire is in the TLL phase, a change of the dot's level causes 
a continuous change in the level occupation, in the total occupation and in 
the free energy. We have shown that this change can be explained within the 
random phase approximation. 

On the other hand, when the wire is described by one of the other couple 
of phases, the FM or the CDW, these thermodynamic quantities behave in 
a different way. There is an abrupt jump in the dot occupation when its 
orbital crosses the chemical potential of the lead. This jump is accompanied 
by an inversion of the occupation of each site in the lead, and in an abrupt 
change of the first derivative of the free energy. We have proven that this is 
a result of a simple level crossing in the FM case, whereas it is a sign of a 
first order quantum phase transition in the case of a CDW. 
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Another influence of tlie dot on tlie coupled system is the Friedel oscil- 
lations (FO) in the wire, which are investigated in chapter HI In metallic 
phases, such as the TLL one, these oscillations decay as a power law, with an 
exponent depending on the interaction strength. If the wire is an insulator, 
as the CDW phase is, the FO decay exponentially. 

Once disorder is introduced in the wire, its phase changes to an Anderson 
insulator (AI). However, for a very weak disorder applied onto a CDW, a 
finite wire may still be described as a Mott insulator (MI). The effect of 
the disorder on the FO decay can be described as another exponential decay 
factor, with a characteristic decay length. We have shown that for a fixed 
weak disorder, an enhancement of the interactions in the AI and the MI 
phases leads to different results: the decay length decreases in the AI phase, 
while it increases in the MI case. This difference was explained according to 
the interplay between interactions and disorder. 

In the AI regime, we have proven that the decay length can be associated 
with Anderson localization length. Our results, presenting a decrease of the 
localization length with increasing interactions, confirm previous predictions. 

The other type of mesoscopic systems we investigate is a 2D QD. Consid- 
ering such a disordered QD, with interacting spinless electrons, we suggest 
in chapter [5] to use the numerical particle-hole DMRG (PH-DMRG) method 
in order to approximate the system's ground state. We have shown that an 
improvement of the PH-DMRG truncation method leads to results which are 
much more accurate than those obtained by Hartree-Fock approximations. 
Furthermore, a significant improvement of the accuracy was exhibited when 
the revised PH-DMRG method was used to calculate the addition spectrum 
of the QD. The suggested method thus opens a door to accurate calculations 
of ground-state properties in two dimensions. 

In chapter [6] we have investigated the lowest states of the 2D QD occu- 
pied by spin 1/2 electrons with interactions and in the presence of spin-orbit 
coupling. We have shown that at certain values of the spin-orbit coupling 
one can obtain a level crossing between the lowest two many-body levels of 
a doubly occupied dot. At the crossing point these two states have identical 
charge distributions, whereas they are different in their spin degree of free- 
dom. Therefore, the level crossing is accompanied by a finite magnetization 
of the ground state, and a finite g-factor is obtained, in contrast to the usual 
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g — case for an even number of electrons. An investigation of the size de- 
pendence of this phenomenon suggests that it might have significant impact 
on g-f actor measurements. 
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